TRestAxionOptics.cxx

/******************** REST disclaimer ***********************************
 * This file is part of the REST software framework.                     *
 *                                                                       *
 * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza)                *
 * For more information see http://gifna.unizar.es/trex                  *
 *                                                                       *
 * REST is free software: you can redistribute it and/or modify          *
 * it under the terms of the GNU General Public License as published by  *
 * the Free Software Foundation, either version 3 of the License, or     *
 * (at your option) any later version.                                   *
 *                                                                       *
 * REST is distributed in the hope that it will be useful,               *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the          *
 * GNU General Public License for more details.                          *
 *                                                                       *
 * You should have a copy of the GNU General Public License along with   *
 * REST in $REST_PATH/LICENSE.                                           *
 * If not, see http://www.gnu.org/licenses/.                             *
 * For the list of contributors see $REST_PATH/CREDITS.                  *
 *************************************************************************/
 
 
#include "TRestAxionOptics.h"
 
#include <TAxis.h>
#include <TGraph.h>
#include <TH1F.h>
#include <TH2F.h>
 
using namespace std;
 
#include "TRestPhysics.h"
using namespace REST_Physics;
 
ClassImp(TRestAxionOptics);
 
TRestAxionOptics::TRestAxionOptics() : TRestMetadata() { Initialize(); }
 
TRestAxionOptics::TRestAxionOptics(const char* cfgFileName, string name) : TRestMetadata(cfgFileName) {
    RESTDebug << "Entering TRestAxionOptics constructor( cfgFileName, name )" << RESTendl;
    RESTDebug << "File: " << cfgFileName << " Name: " << name << RESTendl;
}
 
TRestAxionOptics::~TRestAxionOptics() {
    if (fEntranceMask) delete fEntranceMask;
    if (fExitMask) delete fExitMask;
    if (fMiddleMask) delete fMiddleMask;
}
 
void TRestAxionOptics::Initialize() {
    SetLibraryVersion(LIBRARY_VERSION);
 
    if (fOpticsFile != "") {
        std::string fullPathFileName = SearchFile(fOpticsFile);
 
        if (!TRestTools::ReadASCIITable(fullPathFileName, fOpticsData, 3)) {
            RESTError << "TRestAxionOptics. Error reading optics file : " << fOpticsFile << RESTendl;
        }
 
        // std::cout << "Reading table" << std::endl;
        // TRestTools::PrintTable(fOpticsData, 6, 9);
    }
 
    if (fEntranceMask != nullptr) {
        delete fEntranceMask;
        fEntranceMask = nullptr;
    }
    fEntranceMask = new TRestCombinedMask();
    fEntranceMask->SetName("Entrance");
    fEntranceMask->SetTitle("Optics entrance mask");
    fEntranceMask->SetVerboseLevel(this->GetVerboseLevel());
 
    if (fExitMask != nullptr) {
        delete fExitMask;
        fExitMask = nullptr;
    }
    fExitMask = new TRestCombinedMask();
    fExitMask->SetName("Exit");
    fExitMask->SetTitle("Optics exit mask");
    fExitMask->SetVerboseLevel(this->GetVerboseLevel());
 
    if (fMiddleMask != nullptr) {
        delete fMiddleMask;
        fMiddleMask = nullptr;
    }
    fMiddleMask = new TRestCombinedMask();
    fMiddleMask->SetName("Middle");
    fMiddleMask->SetTitle("Optics middle mask");
    fMiddleMask->SetVerboseLevel(this->GetVerboseLevel());
}
 
Int_t TRestAxionOptics::GetNumberOfReflections() {
    if (fFirstInteraction && fSecondInteraction) return 2;
    if (fFirstInteraction) return 1;
    if (fSecondInteraction) return 1;
    return 0;
}
 
Int_t TRestAxionOptics::TransportToEntrance(const TVector3& pos, const TVector3& dir) {
    RESTDebug << "TRestAxionOptics::TransportToEntrance" << RESTendl;
    if (pos.Z() > GetEntrancePositionZ()) {
        RESTWarning << "TRestAxionOptics::TransportToEntrance" << RESTendl;
        RESTWarning << "The particle should be placed before the entrance!" << RESTendl;
        return 0;
    }
    if (dir.Z() <= 0) {
        RESTWarning << "TRestAxionOptics::TransportToEntrance" << RESTendl;
        RESTWarning << "Photon is not moving on the positive Z-direction!" << RESTendl;
        return 0;
    }
    fEntrancePosition =
        REST_Physics::MoveToPlane(pos, dir, TVector3(0, 0, 1), TVector3(0, 0, GetEntrancePositionZ()));
 
    Double_t x = fEntrancePosition.X();
    Double_t y = fEntrancePosition.Y();
    return fEntranceMask->GetRegion(x, y);
}
 
Int_t TRestAxionOptics::TransportToMiddle(const TVector3& pos, const TVector3& dir) {
    RESTDebug << "TRestAxionOptics::TransportToMiddle" << RESTendl;
    if (pos.Z() > 0 || pos.Z() < GetEntrancePositionZ()) {
        RESTWarning << "TRestAxionOptics::TransportToMiddle" << RESTendl;
        RESTWarning << "The particle should be placed between entrance and middle!" << RESTendl;
        return 0;
    }
    if (dir.Z() <= 0) {
        RESTWarning << "TRestAxionOptics::TransportToMiddle" << RESTendl;
        RESTWarning << "Photon is not moving on the positive Z-direction!" << RESTendl;
        RESTWarning << "Direction : ( " << dir.X() << ", " << dir.Y() << ", " << dir.Z() << ")" << RESTendl;
        return 0;
    }
 
    fMiddlePosition = REST_Physics::MoveToPlane(pos, dir, TVector3(0, 0, 1), TVector3(0, 0, 0));
    fMiddleDirection = dir;
 
    Double_t x = fMiddlePosition.X();
    Double_t y = fMiddlePosition.Y();
    return fMiddleMask->GetRegion(x, y);
}
 
Int_t TRestAxionOptics::TransportToExit(const TVector3& pos, const TVector3& dir) {
    RESTDebug << "TRestAxionOptics::TransportToExit" << RESTendl;
    if (pos.Z() < 0 || pos.Z() > GetExitPositionZ()) {
        RESTWarning << "TRestAxionOptics::TransportToExit" << RESTendl;
        RESTWarning << "The particle should be placed between middle and exit!" << RESTendl;
        return 0;
    }
    if (dir.Z() <= 0) {
        RESTWarning << "TRestAxionOptics::TransportToExit" << RESTendl;
        RESTWarning << "Photon is not moving on the positive Z-direction!" << RESTendl;
        return 0;
    }
 
    fExitPosition =
        REST_Physics::MoveToPlane(pos, dir, TVector3(0, 0, 1), TVector3(0, 0, GetExitPositionZ()));
    fExitDirection = dir;
 
    Double_t x = fExitPosition.X();
    Double_t y = fExitPosition.Y();
    return fExitMask->GetRegion(x, y);
}
 
void TRestAxionOptics::ResetPositions() {
    fOriginPosition = TVector3(0, 0, 0);
 
    fEntrancePosition = TVector3(0, 0, 0);
    fEntranceDirection = TVector3(0, 0, 0);
 
    fFirstInteractionPosition = TVector3(0, 0, 0);
 
    fMiddlePosition = TVector3(0, 0, 0);
    fMiddleDirection = TVector3(0, 0, 0);
 
    fSecondInteractionPosition = TVector3(0, 0, 0);
 
    fExitPosition = TVector3(0, 0, 0);
    fExitDirection = TVector3(0, 0, 0);
 
    fCurrentMirror = -1;
}
 
TVector3 TRestAxionOptics::GetLastGoodPosition() {
    if (GetExitDirection().Mag() > 0)
        return GetExitPosition();
    else if (GetMiddleDirection().Mag() > 0)
        return GetMiddlePosition();
 
    return GetEntrancePosition();
}
 
TVector3 TRestAxionOptics::GetLastGoodDirection() {
    if (GetExitDirection().Mag() > 0)
        return GetExitDirection();
    else if (GetMiddleDirection().Mag() > 0)
        return GetMiddleDirection();
 
    return GetEntranceDirection();
}
 
Double_t TRestAxionOptics::PropagatePhoton(const TVector3& pos, const TVector3& dir, Double_t energy) {
    RESTDebug << " --> Entering TRestAxionOptics::PropagatePhoton" << RESTendl;
 
    RESTDebug << "Reseting positions" << RESTendl;
    ResetPositions();
 
    fOriginPosition = pos;
    fEntranceDirection = dir;
 
    Int_t entranceRegion = TransportToEntrance(fOriginPosition, fEntranceDirection);
    RESTDebug << "Entrance region : " << entranceRegion << RESTendl;
 
    if (entranceRegion == 0) return 0.0;
 
    SetMirror();
 
    FirstMirrorReflection(fEntrancePosition, fEntranceDirection);
 
    Int_t middleRegion = TransportToMiddle(fFirstInteractionPosition, fMiddleDirection);
    RESTDebug << "Middle region : " << middleRegion << RESTendl;
    if (middleRegion != entranceRegion) return 0.0;
 
    SecondMirrorReflection(fMiddlePosition, fMiddleDirection);
 
 
    Int_t exitRegion = TransportToExit(fSecondInteractionPosition, fExitDirection);
    RESTDebug << "Exit region : " << exitRegion << RESTendl;
    if (exitRegion != middleRegion) return 0.0;
 
    return GetPhotonReflectivity(energy);
}
 
void TRestAxionOptics::InitFromConfigFile() {
    if (fMirrorProperties) {
        delete fMirrorProperties;
        fMirrorProperties = nullptr;
    }
 
    fMirrorProperties = (TRestAxionOpticsMirror*)this->InstantiateChildMetadata("TRestAxionOpticsMirror");
 
    TRestMetadata::InitFromConfigFile();
 
    // If we recover the metadata class from a ROOT file we will need to call Initialize ourselves
    this->Initialize();
}
 
void TRestAxionOptics::PrintMetadata() {
    TRestMetadata::PrintMetadata();
 
    RESTMetadata << " - Optics file : " << fOpticsFile << RESTendl;
    RESTMetadata << "---------" << RESTendl;
    RESTMetadata << "Entrance position in Z : " << GetEntrancePositionZ() << " mm" << RESTendl;
    RESTMetadata << "Exit position in Z : " << GetExitPositionZ() << " mm" << RESTendl;
    RESTMetadata << "Low radial limit : " << GetRadialLimits().first << " mm" << RESTendl;
    RESTMetadata << "High radial limit : " << GetRadialLimits().second << " mm" << RESTendl;
    RESTMetadata << "---------" << RESTendl;
    RESTMetadata << " " << RESTendl;
    RESTMetadata << " Use \"this->PrintMasks()\" to get masks info" << RESTendl;
    RESTMetadata << " Use \"this->PrintMirror()\" to get mirror info" << RESTendl;
    RESTMetadata << "+++++++++++++++++++++++++++++++++++++++++++++++++" << RESTendl;
}
 
void TRestAxionOptics::PrintMasks() {
    if (fEntranceMask) fEntranceMask->PrintMetadata();
    if (fMiddleMask) fMiddleMask->PrintMetadata();
    if (fExitMask) fExitMask->PrintMetadata();
}
 
void TRestAxionOptics::PrintMirror() {
    if (fMirrorProperties) fMirrorProperties->PrintMetadata();
}
 
void TRestAxionOptics::PrintEntranceMask() {
    if (fEntranceMask)
        fEntranceMask->PrintMetadata();
    else
        RESTWarning << "TRestAxionOptics::PrintEntranceMask. Not available" << RESTendl;
}
 
void TRestAxionOptics::PrintMiddleMask() {
    if (fMiddleMask)
        fMiddleMask->PrintMetadata();
    else
        RESTWarning << "TRestAxionOptics::PrintMiddleMask. Not available" << RESTendl;
}
 
void TRestAxionOptics::PrintExitMask() {
    if (fExitMask)
        fExitMask->PrintMetadata();
    else
        RESTWarning << "TRestAxionOptics::PrintExitMask. Not available" << RESTendl;
}
 
void TRestAxionOptics::PrintPhotonTrackingSummary() {
    std::cout << "Photon tracking summary" << std::endl;
    std::cout << "-----------------------" << std::endl;
    std::cout << "Origin : ( " << fOriginPosition.X() << ", " << fOriginPosition.Y() << ", "
              << fOriginPosition.Z() << ")" << std::endl;
 
    std::cout << "Entrance Direction : ( " << fEntranceDirection.X() << ", " << fEntranceDirection.Y() << ", "
              << fEntranceDirection.Z() << ")" << std::endl;
 
    std::cout << "Entrance : ( " << fEntrancePosition.X() << ", " << fEntrancePosition.Y() << ", "
              << fEntrancePosition.Z() << ")" << std::endl;
 
    std::cout << "Entrance radius : "
              << TMath::Sqrt(fEntrancePosition.X() * fEntrancePosition.X() +
                             fEntrancePosition.Y() * fEntrancePosition.Y())
              << std::endl;
 
    std::cout << "FirstInteraction : ( " << fFirstInteractionPosition.X() << ", "
              << fFirstInteractionPosition.Y() << ", " << fFirstInteractionPosition.Z() << ")" << std::endl;
 
    std::cout << "Middle Direction : ( " << fMiddleDirection.X() << ", " << fMiddleDirection.Y() << ", "
              << fMiddleDirection.Z() << ")" << std::endl;
 
    std::cout << "Middle : ( " << fMiddlePosition.X() << ", " << fMiddlePosition.Y() << ", "
              << fMiddlePosition.Z() << ")" << std::endl;
 
    std::cout << "Middle radius : "
              << TMath::Sqrt(fMiddlePosition.X() * fMiddlePosition.X() +
                             fMiddlePosition.Y() * fMiddlePosition.Y())
              << std::endl;
 
    std::cout << "SecondInteraction : ( " << fSecondInteractionPosition.X() << ", "
              << fSecondInteractionPosition.Y() << ", " << fSecondInteractionPosition.Z() << ")" << std::endl;
 
    std::cout << "Exit Direction : ( " << fExitDirection.X() << ", " << fExitDirection.Y() << ", "
              << fExitDirection.Z() << ")" << std::endl;
 
    std::cout << "Exit : ( " << fExitPosition.X() << ", " << fExitPosition.Y() << ", " << fExitPosition.Z()
              << ")" << std::endl;
 
    std::cout << "Exit radius : "
              << TMath::Sqrt(fExitPosition.X() * fExitPosition.X() + fExitPosition.Y() * fExitPosition.Y())
              << std::endl;
    std::cout << "-----------------------" << std::endl;
}
 
TPad* TRestAxionOptics::CreatePad(Int_t nx, Int_t ny) {
    if (fPad != nullptr) {
        delete fPad;
        fPad = nullptr;
    }
 
    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);
 
    return fPad;
}
 
Double_t TRestAxionOptics::GetPhotonReflectivity(Double_t energy) {
    if (!fMirrorProperties) return 0;
 
    Double_t reflectivity = 1.;
    if (IsFirstMirrorReflection()) {
        Double_t angle = REST_Physics::GetVectorsAngle(fEntranceDirection, fMiddleDirection);
 
        angle = angle * units("deg") / units("rad") / 2.;  // Incidence angle in degrees. We divide by 2
 
        reflectivity *= fMirrorProperties->GetReflectivity(angle, energy);
    }
    if (IsSecondMirrorReflection()) {
        Double_t angle = REST_Physics::GetVectorsAngle(fMiddleDirection, fExitDirection);
 
        angle = angle * units("deg") / units("rad") / 2.;  // Incidence angle in degrees. We divide by 2
 
        reflectivity *= fMirrorProperties->GetReflectivity(angle, energy);
    }
    return reflectivity;
}
 
TPad* TRestAxionOptics::DrawParticleTracks(Double_t deviation, Int_t particles) {
    DrawMirrors();
 
    for (int n = 0; n < particles; n++) {
        Double_t r = fRandom->Uniform(GetRadialLimits().first, GetRadialLimits().second);
        Double_t angle = fRandom->Uniform(0, 2 * TMath::Pi());
        TVector3 origin(r * TMath::Cos(angle), r * TMath::Sin(angle), -3 * fMirrorLength);
        TVector3 direction(fRandom->Uniform(-deviation, deviation), fRandom->Uniform(-deviation, deviation),
                           1);
        direction = direction.Unit();
 
        Double_t reflectivity = PropagatePhoton(origin, direction, 1);
 
        PrintPhotonTrackingSummary();
 
        TGraph* gr = new TGraph();
 
        Double_t rO = TMath::Sqrt(origin.X() * origin.X() + origin.Y() * origin.Y());
        gr->SetPoint(gr->GetN(), origin.Z(), rO);
 
 
        Double_t rEntrance = TMath::Sqrt(fEntrancePosition.X() * fEntrancePosition.X() +
                                         fEntrancePosition.Y() * fEntrancePosition.Y());
        if (rEntrance > 0) gr->SetPoint(gr->GetN(), fEntrancePosition.Z(), rEntrance);
 
 
        Double_t rFirst = TMath::Sqrt(fFirstInteractionPosition.X() * fFirstInteractionPosition.X() +
                                      fFirstInteractionPosition.Y() * fFirstInteractionPosition.Y());
 
        if (rFirst > 0) gr->SetPoint(gr->GetN(), fFirstInteractionPosition.Z(), rFirst);
 
 
        Double_t rMiddle = TMath::Sqrt(fMiddlePosition.X() * fMiddlePosition.X() +
                                       fMiddlePosition.Y() * fMiddlePosition.Y());
        if (rMiddle > 0) gr->SetPoint(gr->GetN(), fMiddlePosition.Z(), rMiddle);
 
 
        Double_t rSecond = TMath::Sqrt(fSecondInteractionPosition.X() * fSecondInteractionPosition.X() +
                                       fSecondInteractionPosition.Y() * fSecondInteractionPosition.Y());
 
        if (rSecond > 0) gr->SetPoint(gr->GetN(), fSecondInteractionPosition.Z(), rSecond);
 
 
        Double_t rExit =
            TMath::Sqrt(fExitPosition.X() * fExitPosition.X() + fExitPosition.Y() * fExitPosition.Y());
        if (rExit > 0) gr->SetPoint(gr->GetN(), fExitPosition.Z(), rExit);
 
 
        if (reflectivity > 0) {
            TVector3 end = REST_Physics::MoveToPlane(fExitPosition, fExitDirection, TVector3(0, 0, 1),
                                                     TVector3(0, 0, 3 * fMirrorLength));
            Double_t rEnd = TMath::Sqrt(end.X() * end.X() + end.Y() * end.Y());
            if (rEnd > 0) gr->SetPoint(gr->GetN(), end.Z(), rEnd);
        }
 
 
        gr->SetLineWidth(1);
        if (GetMirror() < 0)
            gr->SetLineColor(kBlack);
        else
            gr->SetLineColor(20 + GetMirror() % 20);
 
        gr->Draw("L");
    }
    return fPad;
}
 
Double_t TRestAxionOptics::FindFocal(Double_t from, Double_t to, Double_t energy, Double_t precision,
                                     Bool_t recalculate, Int_t particles) {
    RESTDebug << "Entering TRestAxionOptics::FindFocal" << RESTendl;
 
    if (fFocal > 0 && recalculate == false) return fFocal;
 
    Double_t focal = (from + to) / 2.;
    Double_t spotSize = CalculateSpotSize(energy, focal);
 
    Double_t step = (to - from) / 10.;
 
    if (step < 0) {
        step = from;
        from = to;
        to = step;
        step = (to - from) / 10.;
    }
 
    RESTInfo << "Calculating focal : " << focal << RESTendl;
    for (Double_t f = from; f < to; f += step) {
        Double_t size = CalculateSpotSize(energy, f);
        if (size < spotSize) {
            spotSize = size;
            focal = f;
        }
    }
 
    if (step > precision) return FindFocal(focal - step, focal + step, energy, precision, recalculate);
 
    fFocal = focal;
    return fFocal;
}
 
Double_t TRestAxionOptics::CalculateSpotSize(Double_t energy, Double_t z, Int_t particles) {
    Double_t sum = 0;
    Int_t nSum = 0;
    Double_t deviation = 0;
    for (int n = 0; n < particles; n++) {
        Double_t reflectivity = PropagateMonteCarloPhoton(energy, deviation);
 
        if (reflectivity == 0) continue;
 
        TVector3 posZ =
            REST_Physics::MoveToPlane(fExitPosition, fExitDirection, TVector3(0, 0, 1), TVector3(0, 0, z));
        Double_t x = posZ.X();
        Double_t y = posZ.Y();
 
        sum += reflectivity * x * x + y * y;
        nSum++;
    }
 
    if (nSum > 0) sum = TMath::Sqrt(sum / nSum);
 
    return sum;
}
 
Double_t TRestAxionOptics::PropagateMonteCarloPhoton(Double_t energy, Double_t deviation) {
    Double_t x = fRandom->Uniform(0, 1.5 * GetRadialLimits().second);
    Double_t y = fRandom->Uniform(0, 1.5 * GetRadialLimits().second);
    Double_t r2 = x * x + y * y;
    while (r2 > 1.5 * GetRadialLimits().second || r2 < 0.5 * GetRadialLimits().second) {
        x = fRandom->Uniform(0, 1.5 * GetRadialLimits().second);
        y = fRandom->Uniform(0, 1.5 * GetRadialLimits().second);
        r2 = x * x + y * y;
    }
 
    Double_t r = fRandom->Uniform(0.5 * GetRadialLimits().first, 1.5 * GetRadialLimits().second);
    Double_t angle = fRandom->Uniform(0, 2 * TMath::Pi());
    TVector3 origin(r * TMath::Cos(angle), r * TMath::Sin(angle), -3 * fMirrorLength);
 
    Double_t theta = fRandom->Uniform(0, deviation);
    Double_t phi = fRandom->Uniform(0, 2 * TMath::Pi());
 
    TVector3 direction(0, 0, 1);
    direction.SetTheta(theta);
    direction.SetPhi(phi);
    direction.SetMag(1.);
 
    return PropagatePhoton(origin, direction, energy);
}
 
TPad* TRestAxionOptics::DrawScatterMaps(Double_t z, Double_t energy, Double_t deviation, Int_t particles,
                                        Double_t focalHint) {
    TRestAxionOptics::CreatePad(2, 2);
 
    fPad->cd();
 
    TGraph* grEntrance = new TGraph();
    TGraph* grExit = new TGraph();
    TGraph* grZ = new TGraph();
    TGraph* grFocal = new TGraph();
 
    Double_t focal = FindFocal(focalHint - 500, focalHint + 500, energy, 1);
 
    for (int n = 0; n < particles; n++) {
        Double_t reflectivity = PropagateMonteCarloPhoton(energy, deviation);
 
        if (fFirstInteractionPosition.Z() == 0) continue;  // The photon hits the entrance mask
        grEntrance->SetPoint(grEntrance->GetN(), fEntrancePosition.X(), fEntrancePosition.Y());
 
        if (reflectivity == 0) continue;  // The photon hits any mask
        grExit->SetPoint(grExit->GetN(), fExitPosition.X(), fExitPosition.Y());
 
        TVector3 posZ =
            REST_Physics::MoveToPlane(fExitPosition, fExitDirection, TVector3(0, 0, 1), TVector3(0, 0, z));
        grZ->SetPoint(grZ->GetN(), posZ.X(), posZ.Y());
 
        TVector3 posFocal = REST_Physics::MoveToPlane(fExitPosition, fExitDirection, TVector3(0, 0, 1),
                                                      TVector3(0, 0, focal));
        grFocal->SetPoint(grFocal->GetN(), posFocal.X(), posFocal.Y());
    }
 
    fPad->cd(1);
 
    grEntrance->GetXaxis()->SetLimits(-GetRadialLimits().second * 1.15, GetRadialLimits().second * 1.15);
    grEntrance->GetHistogram()->SetMaximum(GetRadialLimits().second * 1.15);
    grEntrance->GetHistogram()->SetMinimum(-GetRadialLimits().second * 1.15);
 
    grEntrance->SetTitle("Entrance plane");
    grEntrance->GetXaxis()->SetTitle("X [mm]");
    grEntrance->GetXaxis()->SetTitleSize(0.04);
    grEntrance->GetXaxis()->SetLabelSize(0.04);
    grEntrance->GetXaxis()->SetNdivisions(5);
    grEntrance->GetYaxis()->SetTitle("Y [mm]");
    grEntrance->GetYaxis()->SetTitleOffset(1.4);
    grEntrance->GetYaxis()->SetTitleSize(0.04);
    grEntrance->GetYaxis()->SetLabelSize(0.04);
 
    grEntrance->SetMarkerStyle(1);
    grEntrance->Draw("AP");
 
    fPad->cd(2);
 
    grExit->SetTitle("Exit plane");
    grExit->GetXaxis()->SetTitle("X [mm]");
    grExit->GetXaxis()->SetTitleSize(0.04);
    grExit->GetXaxis()->SetLabelSize(0.04);
    grExit->GetXaxis()->SetNdivisions(5);
    grExit->GetYaxis()->SetTitle("Y [mm]");
    grExit->GetYaxis()->SetTitleOffset(1.4);
    grExit->GetYaxis()->SetTitleSize(0.04);
    grExit->GetYaxis()->SetLabelSize(0.04);
 
    grExit->SetMarkerStyle(1);
    grExit->Draw("AP");
 
    fPad->cd(3);
 
    Double_t xMin = TMath::MinElement(grZ->GetN(), grZ->GetX());
    Double_t xMax = TMath::MaxElement(grZ->GetN(), grZ->GetX());
    Double_t yMin = TMath::MinElement(grZ->GetN(), grZ->GetY());
    Double_t yMax = TMath::MaxElement(grZ->GetN(), grZ->GetY());
 
    std::string zTitle = "User plane. Z = " + DoubleToString(z) + "mm";
    grZ->SetTitle(zTitle.c_str());
    grZ->GetXaxis()->SetLimits(1.5 * xMin, 1.5 * xMax);
    grZ->GetHistogram()->SetMaximum(1.5 * yMax);
    grZ->GetHistogram()->SetMinimum(1.5 * yMin);
 
    grZ->GetXaxis()->SetTitle("X [mm]");
    grZ->GetXaxis()->SetTitleSize(0.04);
    grZ->GetXaxis()->SetLabelSize(0.04);
    grZ->GetXaxis()->SetNdivisions(5);
    grZ->GetYaxis()->SetTitle("Y [mm]");
    grZ->GetYaxis()->SetTitleOffset(1.4);
    grZ->GetYaxis()->SetTitleSize(0.04);
    grZ->GetYaxis()->SetLabelSize(0.04);
 
    grZ->SetMarkerStyle(1);
    grZ->Draw("AP");
 
    fPad->cd(4);
 
    std::string focalTitle = "Focal plane. Z = " + DoubleToString(focal) + "mm";
    grFocal->SetTitle(focalTitle.c_str());
    grFocal->GetXaxis()->SetLimits(-10, 10);
    grFocal->GetHistogram()->SetMaximum(10);
    grFocal->GetHistogram()->SetMinimum(-10);
 
    grFocal->GetXaxis()->SetTitle("X [mm]");
    grFocal->GetXaxis()->SetTitleSize(0.04);
    grFocal->GetXaxis()->SetLabelSize(0.04);
    grFocal->GetXaxis()->SetNdivisions(5);
    grFocal->GetYaxis()->SetTitle("Y [mm]");
    grFocal->GetYaxis()->SetTitleOffset(1.4);
    grFocal->GetYaxis()->SetTitleSize(0.04);
    grFocal->GetYaxis()->SetLabelSize(0.04);
 
    grFocal->SetMarkerStyle(1);
    grFocal->Draw("AP");
 
    return fPad;
}
 
TPad* TRestAxionOptics::DrawDensityMaps(Double_t z, Double_t energy, Double_t deviation, Int_t particles,
                                        Double_t focalHint) {
    Double_t focal = FindFocal(focalHint - 500, focalHint + 500, energy, 1);
 
    TRestAxionOptics::CreatePad(2, 2);
 
    fPad->cd();
 
    Double_t lowL = -1.15 * GetRadialLimits().second;
    Double_t highL = 1.15 * GetRadialLimits().second;
    TH2F* hEntrance = new TH2F("entranceH", "Entrance plane", 500, lowL, highL, 500, lowL, highL);
    TH2F* hExit = new TH2F("exitH", "Exit plane", 500, lowL, highL, 500, lowL, highL);
 
    std::string zTitle = "User plane. Z = " + DoubleToString(z) + "mm";
    TH2F* hZ = new TH2F("zH", zTitle.c_str(), 500, lowL, highL, 500, lowL, highL);
 
    std::string focalTitle = "Focal plane. Z = " + DoubleToString(focal) + "mm";
    TH2F* hFocal = new TH2F("focalH", focalTitle.c_str(), 500, -10, 10, 500, -10, 10);
 
    for (int n = 0; n < particles; n++) {
        Double_t reflectivity = PropagateMonteCarloPhoton(energy, deviation);
 
        if (fFirstInteractionPosition.Z() == 0) continue;  // The photon hits the entrance mask
        hEntrance->Fill(fEntrancePosition.X(), fEntrancePosition.Y());
 
        if (reflectivity == 0) continue;  // The photon hits any mask
        hExit->Fill(fExitPosition.X(), fExitPosition.Y());
 
        TVector3 posZ =
            REST_Physics::MoveToPlane(fExitPosition, fExitDirection, TVector3(0, 0, 1), TVector3(0, 0, z));
        hZ->Fill(posZ.X(), posZ.Y());
 
        TVector3 posFocal = REST_Physics::MoveToPlane(fExitPosition, fExitDirection, TVector3(0, 0, 1),
                                                      TVector3(0, 0, focal));
        hFocal->Fill(posFocal.X(), posFocal.Y());
    }
 
    fPad->cd(1);
 
    hEntrance->GetXaxis()->SetTitle("X [mm]");
    hEntrance->GetXaxis()->SetTitleSize(0.04);
    hEntrance->GetXaxis()->SetLabelSize(0.04);
    hEntrance->GetXaxis()->SetNdivisions(5);
    hEntrance->GetYaxis()->SetTitle("Y [mm]");
    hEntrance->GetYaxis()->SetTitleOffset(1.4);
    hEntrance->GetYaxis()->SetTitleSize(0.04);
    hEntrance->GetYaxis()->SetLabelSize(0.04);
 
    hEntrance->Draw("colz");
 
    fPad->cd(2);
 
    hExit->GetXaxis()->SetTitle("X [mm]");
    hExit->GetXaxis()->SetTitleSize(0.04);
    hExit->GetXaxis()->SetLabelSize(0.04);
    hExit->GetXaxis()->SetNdivisions(5);
    hExit->GetYaxis()->SetTitle("Y [mm]");
    hExit->GetYaxis()->SetTitleOffset(1.4);
    hExit->GetYaxis()->SetTitleSize(0.04);
    hExit->GetYaxis()->SetLabelSize(0.04);
 
    hExit->Draw("colz");
 
    fPad->cd(3);
 
    hZ->GetXaxis()->SetTitle("X [mm]");
    hZ->GetXaxis()->SetTitleSize(0.04);
    hZ->GetXaxis()->SetLabelSize(0.04);
    hZ->GetXaxis()->SetNdivisions(5);
    hZ->GetYaxis()->SetTitle("Y [mm]");
    hZ->GetYaxis()->SetTitleOffset(1.4);
    hZ->GetYaxis()->SetTitleSize(0.04);
    hZ->GetYaxis()->SetLabelSize(0.04);
 
    hZ->Draw("colz");
 
    fPad->cd(4);
 
    hFocal->GetXaxis()->SetTitle("X [mm]");
    hFocal->GetXaxis()->SetTitleSize(0.04);
    hFocal->GetXaxis()->SetLabelSize(0.04);
    hFocal->GetXaxis()->SetNdivisions(5);
    hFocal->GetYaxis()->SetTitle("Y [mm]");
    hFocal->GetYaxis()->SetTitleOffset(1.4);
    hFocal->GetYaxis()->SetTitleSize(0.04);
    hFocal->GetYaxis()->SetLabelSize(0.04);
 
    hFocal->Draw("colz");
 
    return fPad;
}