Skip to content

Commit

Permalink
[pre-commit.ci] auto fixes from pre-commit.com hooks
Browse files Browse the repository at this point in the history 
for more information, see https://pre-commit.ci
  • Loading branch information
@pre-commit-ci
pre-commit-ci[bot] committed May 3, 2024
1 parent 55bb571 commit 74793f3
Show file tree
Hide file tree
Showing 3 changed files with 66 additions and 74 deletions.
4 changes: 2 additions & 2 deletions inc/TRestAxionMagneticField.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -161,8 +161,8 @@ class TRestAxionMagneticField : public TRestMetadata {
std::vector<Double_t> GetTransversalComponentAlongPath(TVector3 from, TVector3 to, Double_t dl = 1.,
Int_t Nmax = 0);

std::vector<Double_t> GetComponentAlongPath( Int_t axis, TVector3 from, TVector3 to, Double_t dl = 1.,
Int_t Nmax = 0);
std::vector<Double_t> GetComponentAlongPath(Int_t axis, TVector3 from, TVector3 to, Double_t dl = 1.,
Int_t Nmax = 0);

Double_t GetTransversalFieldAverage(TVector3 from, TVector3 to, Double_t dl = 1., Int_t Nmax = 0);

Expand Down
34 changes: 15 additions & 19 deletions src/TRestAxionField.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -210,7 +210,6 @@

#include <TComplex.h>
#include <TVectorD.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>

Expand Down Expand Up @@ -484,7 +483,7 @@ std::pair<Double_t, Double_t> TRestAxionField::GammaTransmissionFieldMapProbabil
Double_t accuracy,
Int_t num_intervals,
Int_t qawo_levels) {
gsl_set_error_handler_off();
gsl_set_error_handler_off();

if (!fMagneticField) {
RESTError << "TRestAxionField::GammaTransmissionFieldMapProbability requires a magnetic field map!"
Expand All @@ -493,8 +492,7 @@ std::pair<Double_t, Double_t> TRestAxionField::GammaTransmissionFieldMapProbabil
return {0.0, 0.0};
}

if( fMagneticField->GetTrackLength() <= 0 )
return {0.0, 0.0};
if (fMagneticField->GetTrackLength() <= 0) return {0.0, 0.0};

double photonMass = 0; // Vacuum
if (fBufferGas) photonMass = fBufferGas->GetPhotonMass(Ea);
Expand Down Expand Up @@ -562,11 +560,10 @@ std::pair<Double_t, Double_t> TRestAxionField::ComputeResonanceIntegral(Double_t

auto start = std::chrono::system_clock::now();

int status = gsl_integration_qag(&F, 0, fMagneticField->GetTrackLength(), accuracy, accuracy, num_intervals,
GSL_INTEG_GAUSS61, workspace, &reprob, &rerr);
int status = gsl_integration_qag(&F, 0, fMagneticField->GetTrackLength(), accuracy, accuracy,
num_intervals, GSL_INTEG_GAUSS61, workspace, &reprob, &rerr);

if( status > 0 )
return {0,status};
if (status > 0) return {0, status};

auto end = std::chrono::system_clock::now();
auto seconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
Expand Down Expand Up @@ -628,20 +625,19 @@ std::pair<Double_t, Double_t> TRestAxionField::ComputeOffResonanceIntegral(Doubl

gsl_integration_qawo_table* wf =
gsl_integration_qawo_table_alloc(q, fMagneticField->GetTrackLength(), GSL_INTEG_COSINE, qawo_levels);
int status = gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &reprob, &rerr);
if( status > 0 )
{
gsl_integration_qawo_table_free(wf);
return {0,status};
}
int status =
gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &reprob, &rerr);
if (status > 0) {
gsl_integration_qawo_table_free(wf);
return {0, status};
}

gsl_integration_qawo_table_set(wf, q, fMagneticField->GetTrackLength(), GSL_INTEG_SINE);
status = gsl_integration_qawo(&F, 0, accuracy, accuracy, num_intervals, workspace, wf, &improb, &imerr);
if( status > 0 )
{
gsl_integration_qawo_table_free(wf);
return {0,status};
}
if (status > 0) {
gsl_integration_qawo_table_free(wf);
return {0, status};
}

auto end = std::chrono::system_clock::now();
auto seconds = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
Expand Down
102 changes: 49 additions & 53 deletions src/TRestAxionMagneticField.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -589,56 +589,56 @@ TCanvas* TRestAxionMagneticField::DrawHistogram(TString projection, TString Bcom
/// are drawm together with their corresponding field intensity profile along the Z-axis.
///
TCanvas* TRestAxionMagneticField::DrawComponents(Int_t volId) {
Int_t divisions = 100;
Int_t divisions = 100;
if (fCanvas != NULL) {
delete fCanvas;
fCanvas = NULL;
}
fCanvas = new TCanvas("fCanvas", "", 1600, 600);

TPad* pad1 = new TPad("pad1", "This is pad1", 0.01, 0.02, 0.99, 0.97);
// pad1->Divide(2, 1);
// pad1->Divide(2, 1);
pad1->Draw();
pad1->cd();

Int_t n = 0;
Int_t n = 0;
Double_t genPositionZ = fPositions[volId][2] - fBoundMax[volId].Z() - 2000;
Double_t finalPositionZ = fPositions[volId][2] + fBoundMax[volId].Z() + 2000;
TVector3 position(0, 0, genPositionZ);
TVector3 direction (0,0,1);

RESTDebug << RESTendl;
RESTDebug << "Start moving along" << RESTendl;
RESTDebug << "++++++++++++++++++" << RESTendl;

TGraph* fieldGr = new TGraph();
Double_t posZ = fPositions[volId][2] - fBoundMax[volId].Z() - 10;
Double_t delta = fBoundMax[volId][2] * 2. / divisions;

while (posZ <= fPositions[volId][2] + fBoundMax[volId].Z()) {
TVector3 posAlongAxis = TVector3(fPositions[volId][0], fPositions[volId][1], posZ);

position = MoveToPlane(position, direction, TVector3(0, 0, 1), posAlongAxis);
Double_t fieldY = this->GetMagneticField(position).Y();

fieldGr->SetPoint(fieldGr->GetN(), posZ, fieldY);

posZ += delta;
}

fieldGr->SetLineWidth(3);
fieldGr->SetLineColor(38 + n);
fieldGr->GetXaxis()->SetLimits(genPositionZ - 500, finalPositionZ + 500);
fieldGr->GetHistogram()->SetMaximum(2.5);
fieldGr->GetHistogram()->SetMinimum(0);
fieldGr->GetXaxis()->SetTitle("Z [mm]");
fieldGr->GetXaxis()->SetTitleSize(0.05);
fieldGr->GetXaxis()->SetLabelSize(0.05);
fieldGr->GetYaxis()->SetTitle("B [T]");
fieldGr->GetYaxis()->SetTitleOffset(1.3);
fieldGr->GetYaxis()->SetTitleSize(0.05);
fieldGr->GetYaxis()->SetLabelSize(0.05);
fieldGr->Draw("AL");
TVector3 position(0, 0, genPositionZ);
TVector3 direction(0, 0, 1);

RESTDebug << RESTendl;
RESTDebug << "Start moving along" << RESTendl;
RESTDebug << "++++++++++++++++++" << RESTendl;

TGraph* fieldGr = new TGraph();
Double_t posZ = fPositions[volId][2] - fBoundMax[volId].Z() - 10;
Double_t delta = fBoundMax[volId][2] * 2. / divisions;

while (posZ <= fPositions[volId][2] + fBoundMax[volId].Z()) {
TVector3 posAlongAxis = TVector3(fPositions[volId][0], fPositions[volId][1], posZ);

position = MoveToPlane(position, direction, TVector3(0, 0, 1), posAlongAxis);
Double_t fieldY = this->GetMagneticField(position).Y();

fieldGr->SetPoint(fieldGr->GetN(), posZ, fieldY);

posZ += delta;
}

fieldGr->SetLineWidth(3);
fieldGr->SetLineColor(38 + n);
fieldGr->GetXaxis()->SetLimits(genPositionZ - 500, finalPositionZ + 500);
fieldGr->GetHistogram()->SetMaximum(2.5);
fieldGr->GetHistogram()->SetMinimum(0);
fieldGr->GetXaxis()->SetTitle("Z [mm]");
fieldGr->GetXaxis()->SetTitleSize(0.05);
fieldGr->GetXaxis()->SetLabelSize(0.05);
fieldGr->GetYaxis()->SetTitle("B [T]");
fieldGr->GetYaxis()->SetTitleOffset(1.3);
fieldGr->GetYaxis()->SetTitleSize(0.05);
fieldGr->GetYaxis()->SetLabelSize(0.05);
fieldGr->Draw("AL");

return fCanvas;
}
Expand Down Expand Up @@ -1286,21 +1286,20 @@ std::vector<Double_t> TRestAxionMagneticField::GetTransversalComponentAlongPath(
/// condition.
///
std::vector<Double_t> TRestAxionMagneticField::GetComponentAlongPath(Int_t axis, TVector3 from, TVector3 to,
Double_t dl, Int_t Nmax) {
Double_t dl, Int_t Nmax) {
std::vector<Double_t> Bt;
if( axis != 0 && axis != 1 && axis != 2 )
{
RESTError << "TRestAxionMagneticField::GetComponentAlongPath. Axis must take values 0,1 or 2" << RESTendl;
return Bt;
}
if (axis != 0 && axis != 1 && axis != 2) {
RESTError << "TRestAxionMagneticField::GetComponentAlongPath. Axis must take values 0,1 or 2"
<< RESTendl;
return Bt;
}
Double_t length = (to - from).Mag();

Double_t diff = dl;
if (Nmax > 0) {
if (length / dl > Nmax) {
diff = length / Nmax;
RESTWarning << "TRestAxionMagneticField::GetComponentAlongPath. Nmax reached!"
<< RESTendl;
RESTWarning << "TRestAxionMagneticField::GetComponentAlongPath. Nmax reached!" << RESTendl;
RESTWarning << "Nmax = " << Nmax << RESTendl;
RESTWarning << "Adjusting differential step to : " << diff << " mm" << RESTendl;
}
Expand All @@ -1309,12 +1308,9 @@ std::vector<Double_t> TRestAxionMagneticField::GetComponentAlongPath(Int_t axis,
TVector3 direction = (to - from).Unit();

for (Double_t d = 0; d < length; d += diff) {
if( axis == 0 )
Bt.push_back(GetMagneticField(from + d * direction).X());
if( axis == 1 )
Bt.push_back(GetMagneticField(from + d * direction).Y());
if( axis == 2 )
Bt.push_back(GetMagneticField(from + d * direction).Z());
if (axis == 0) Bt.push_back(GetMagneticField(from + d * direction).X());
if (axis == 1) Bt.push_back(GetMagneticField(from + d * direction).Y());
if (axis == 2) Bt.push_back(GetMagneticField(from + d * direction).Z());
}

return Bt;
Expand All @@ -1331,7 +1327,7 @@ void TRestAxionMagneticField::SetTrack(const TVector3& position, const TVector3&
fTrackStart = TVector3(0, 0, 0);
fTrackDirection = TVector3(0, 0, 0);
fTrackLength = 0;
return;
return;
}

fTrackStart = trackBounds[0];
Expand Down

0 comments on commit 74793f3

Please sign in to comment.