TRestAxionWolterOptics.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 "TRestAxionWolterOptics.h"
 
using namespace std;
#include <TAxis.h>
#include <TGraph.h>
#include <TH1F.h>
 
#include <cmath>
 
#include "TRestPhysics.h"
using namespace REST_Physics;
ClassImp(TRestAxionWolterOptics);
 
TRestAxionWolterOptics::TRestAxionWolterOptics() : TRestAxionOptics() { Initialize(); }
 
TRestAxionWolterOptics::~TRestAxionWolterOptics() {}
 
TRestAxionWolterOptics::TRestAxionWolterOptics(const char* cfgFileName, string name)
    : TRestAxionOptics(cfgFileName) {
    RESTDebug << "Entering TRestAxionWolterOptics constructor( cfgFileName, name )" << RESTendl;
 
    Initialize();
 
    LoadConfigFromFile(fConfigFileName, name);
 
    if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata();
}
 
void TRestAxionWolterOptics::Initialize() {
    TRestAxionOptics::Initialize();
 
    SetSectionName(this->ClassName());
    SetLibraryVersion(LIBRARY_VERSION);
 
    fR1 = GetR1();
    fR2 = GetR2();
    fR3 = GetR3();
    fR4 = GetR4();
    fR5 = GetR5();
    fAlpha = GetAlpha();
    fThickness = GetThickness();
 
    fXSep.clear();
    for (unsigned int n = 0; n < fAlpha.size(); n++)
        if (fR2[n] == fR3[n]) {
            fXSep.push_back(0);
        } else {
            fXSep.push_back(2 * (fR1[n] - fR3[n] - fMirrorLength * TMath::Sin(fAlpha[n])) /
                            TMath::Tan(fAlpha[n]));
        }
 
    if (fAlpha.size() == 0) return;
 
    fCosAlpha.clear();
    for (const auto& a : fAlpha) fCosAlpha.push_back(TMath::Cos(a));
 
    fCosAlpha_3.clear();
    for (const auto& a : fAlpha) fCosAlpha_3.push_back(TMath::Cos(3 * a));
 
    fFrontVertex.clear();
    for (unsigned int n = 0; n < fAlpha.size(); n++) fFrontVertex.push_back(fR3[n] / TMath::Tan(fAlpha[n]));
 
    fBackVertex.clear();
    for (unsigned int n = 0; n < fAlpha.size(); n++)
        fBackVertex.push_back(fR3[n] / TMath::Tan(3. * fAlpha[n]));
 
    if (fEntranceRingsMask) {
        delete fEntranceRingsMask;
        fEntranceRingsMask = nullptr;
    }
    fEntranceRingsMask = new TRestRingsMask();
 
    std::vector<Double_t> inner, outer;
    for (unsigned int n = 0; n < fR1.size() - 1; n++) {
        inner.push_back(fR1[n] + fThickness[n]);
        outer.push_back(fR1[n + 1]);
    }
    fEntranceRingsMask->SetRadii(inner, outer);
 
    fEntranceMask->AddMask(fSpiderMask);
    fEntranceMask->AddMask(fEntranceRingsMask);
 
    if (fMiddleRingsMask) {
        delete fMiddleRingsMask;
        fMiddleRingsMask = nullptr;
    }
    fMiddleRingsMask = new TRestRingsMask();
 
    inner.clear();
    outer.clear();
    for (unsigned int n = 0; n < fR3.size() - 1; n++) {
        inner.push_back(fR3[n] + fThickness[n]);
        outer.push_back(fR3[n + 1]);
    }
    fMiddleRingsMask->SetRadii(inner, outer);
 
    fMiddleMask->AddMask(fSpiderMask);
    fMiddleMask->AddMask(fMiddleRingsMask);
 
    if (fExitRingsMask) {
        delete fExitRingsMask;
        fExitRingsMask = nullptr;
    }
    fExitRingsMask = new TRestRingsMask();
 
    inner.clear();
    outer.clear();
    for (unsigned int n = 0; n < fR5.size() - 1; n++) {
        inner.push_back(fR5[n] + fThickness[n]);
        outer.push_back(fR5[n + 1]);
    }
    fExitRingsMask->SetRadii(inner, outer);
 
    fExitMask->AddMask(fSpiderMask);
    fExitMask->AddMask(fExitRingsMask);
 
    if (fRandom != nullptr) {
        delete fRandom;
        fRandom = nullptr;
    }
    fRandom = new TRandom3(0);
}
 
Int_t TRestAxionWolterOptics::FirstMirrorReflection(const TVector3& pos, const TVector3& dir) {
    Int_t mirror = GetMirror();
    RESTDebug << "--> Entering TRestAxionWolterOptics::FirstMirrorReflection" << RESTendl;
    RESTDebug << "Mirror: " << mirror << RESTendl;
    if (mirror < 0) {
        RESTError << "TRestAxionWolterOptics::FirstMirrorReflection. Mirror index cannot be negative!"
                  << RESTendl;
        return -1;
    }
 
    if (mirror >= 0 && (unsigned int)mirror >= fFrontVertex.size()) {
        RESTError << "TRestAxionWolterOptics::FirstMirrorReflection. Mirror index above number of mirrors!"
                  << RESTendl;
        return -1;
    }
 
    RESTDebug << "Vertex Z: " << fFrontVertex[mirror] << RESTendl;
    RESTDebug << "Cos Alpha: " << fCosAlpha[mirror] << RESTendl;
    TVector3 vertex(0, 0, fFrontVertex[mirror]);
    Double_t cosA = fCosAlpha[mirror];
 
    fFirstInteractionPosition =
        pos + dir * REST_Physics::GetConeVectorIntersection(pos, dir, TVector3(0, 0, -1), vertex, cosA);
 
    if (fFirstInteractionPosition.Z() < GetEntrancePositionZ() ||
        fFirstInteractionPosition.Z() > -(0.5 * fXSep[mirror])) {
        RESTDebug << "TRestAxionWolterOptics::FirstMirrorReflection. No interaction!" << RESTendl;
        fFirstInteractionPosition = REST_Physics::MoveByDistance(pos, dir, fMirrorLength / 2.);
        fMiddleDirection = fEntranceDirection;
        fFirstInteraction = false;
        return 0;
    }
 
    TVector3 coneNormal = REST_Physics::GetConeNormal(fFirstInteractionPosition, fAlpha[mirror]);
    RESTDebug << "Cone normal: (" << coneNormal.X() << ", " << coneNormal.Y() << ", " << coneNormal.Z() << ")"
              << RESTendl;
 
    fMiddleDirection = GetVectorReflection(fEntranceDirection, coneNormal);
 
    RESTDebug << "<-- Exiting TRestAxionWolterOptics::FirstMirrorReflection" << RESTendl;
 
    fFirstInteraction = true;
    return 1;
}
 
Int_t TRestAxionWolterOptics::SecondMirrorReflection(const TVector3& pos, const TVector3& dir) {
    Int_t mirror = GetMirror();
    if (mirror < 0) {
        RESTError << "TRestAxionWolterOptics::FirstMirrorReflection. Mirror index cannot be negative!"
                  << RESTendl;
        return 0;
    }
 
    if (mirror >= 0 && (unsigned int)mirror >= fFrontVertex.size()) {
        RESTError << "TRestAxionWolterOptics::FirstMirrorReflection. Mirror index above number of mirrors!"
                  << RESTendl;
        return 0;
    }
 
    TVector3 vertex(0, 0, fBackVertex[mirror]);
    Double_t cosA = fCosAlpha_3[mirror];
 
    fSecondInteractionPosition =
        pos + dir * REST_Physics::GetConeVectorIntersection(pos, dir, TVector3(0, 0, -1), vertex, cosA);
 
    if (fSecondInteractionPosition.Z() > GetExitPositionZ() ||
        fSecondInteractionPosition.Z() < (0.5 * fXSep[mirror])) {
        RESTDebug << "TRestAxionWolterOptics::SecondMirrorReflection. No interaction!" << RESTendl;
        fSecondInteractionPosition = REST_Physics::MoveByDistance(pos, dir, fMirrorLength / 2.);
        fExitDirection = fMiddleDirection;
        fSecondInteraction = false;
        return 0;
    }
 
    TVector3 coneNormal = REST_Physics::GetConeNormal(fSecondInteractionPosition, 3 * fAlpha[mirror]);
 
    fExitDirection = GetVectorReflection(fMiddleDirection, coneNormal);
 
    fSecondInteraction = true;
    return 1;
}
 
void TRestAxionWolterOptics::InitFromConfigFile() {
    if (fSpiderMask) {
        delete fSpiderMask;
        fSpiderMask = nullptr;
    }
 
    fSpiderMask = (TRestSpiderMask*)this->InstantiateChildMetadata("TRestSpiderMask");
 
    if (fSpiderMask == nullptr) {
        RESTWarning << "TRestAxionWolterOptics requires usually a TRestSpiderMask definition" << RESTendl;
    } else {
    }
 
    TRestAxionOptics::InitFromConfigFile();
 
    // If we recover the metadata class from ROOT file we will need to call Initialize ourselves
    this->Initialize();
}
 
void TRestAxionWolterOptics::PrintParameters() {
    if (fR3.size() > 0) {
        for (unsigned int n = 0; n < fR3.size(); n++) {
            Double_t dR1 = fR1[n] - fR3[n] - fMirrorLength * TMath::Sin(fAlpha[n]);
            Double_t dR5 = fR5[n] - fR3[n] + fMirrorLength * TMath::Sin(3 * fAlpha[n]);
            if (n % 10 == 0)
                std::cout << "## R1\tdelta R1\tR3\tR5\tdelta R5\talpha\tCosAlpha\tFrontVertex\tBackVertex"
                          << std::endl;
 
            std::cout << fR1[n] << "\t" << dR1 << "\t" << fR3[n] << "\t" << fR5[n] << "\t" << dR5 << "\t"
                      << fAlpha[n] << "\t" << fCosAlpha[n] << "\t" << fFrontVertex[n] << "\t"
                      << fBackVertex[n] << std::endl;
        }
    }
}
 
void TRestAxionWolterOptics::PrintSpider() {
    if (fSpiderMask) fSpiderMask->PrintMetadata();
}
 
void TRestAxionWolterOptics::PrintMetadata() { TRestAxionOptics::PrintMetadata(); }
 
TPad* TRestAxionWolterOptics::DrawMirrors() {
    TRestAxionOptics::CreatePad();
 
    fPad->cd();
    std::vector<TGraph*> graphCollection;
    for (unsigned int mirror = 0; mirror < fR3.size(); mirror++) {
        TGraph* gr = new TGraph();  //"Mirror" + IntegerToString(mirror + 1));
 
        Double_t lX = fMirrorLength * fCosAlpha[mirror];
 
        gr->SetPoint(0, -lX, fR1[mirror]);
        gr->SetPoint(1, 0, fR3[mirror]);
        gr->SetPoint(2, lX, fR5[mirror]);
 
        gr->GetXaxis()->SetLimits(-3.5 * lX, 3.5 * lX);
        gr->GetHistogram()->SetMaximum(fR1.back() * 1.15);
        gr->GetHistogram()->SetMinimum(fR1.front() * 0.8);
 
        gr->GetXaxis()->SetTitle("Z [mm]");
        gr->GetXaxis()->SetTitleSize(0.04);
        gr->GetXaxis()->SetLabelSize(0.04);
        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->SetLineWidth(6 * fThickness[mirror]);
        gr->SetLineColor(20 + mirror % 20);
        if (mirror == 0)
            gr->Draw("AL");
        else
            gr->Draw("L");
    }
 
    return fPad;
}