TRestAxionMagneticField.h

/*************************************************************************
 * 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.                  *
 *************************************************************************/

#ifndef _TRestAxionMagneticField
#define _TRestAxionMagneticField

#include <TRestMetadata.h>
#include <iostream>

#include "TCanvas.h"
#include "TH2D.h"
#include "TVector3.h"
#include "TVectorD.h"

#include "TRestMesh.h"

struct MagneticFieldVolume {
    TVector3 position;

    TRestMesh mesh;

    std::vector<std::vector<std::vector<TVector3>>> field;
};

class TRestAxionMagneticField : public TRestMetadata {
   private:
    std::vector<std::string> fFileNames;  //<

    std::vector<TVector3> fPositions;  //<

    std::vector<TVector3> fConstantField;  //<

    std::vector<TVector3> fMeshSize;  //<

    std::vector<TString> fMeshType;  //<

    std::vector<TVector3> fBoundMax;  //<

    std::vector<MagneticFieldVolume> fMagneticFieldVolumes;  

    TH2D* fHisto;  

    TCanvas* fCanvas;  

    void Initialize();

    void InitFromConfigFile();

    void LoadMagneticFieldData(MagneticFieldVolume& mVol, std::vector<std::vector<Float_t>> data);

    TVector3 GetMagneticVolumeNode(MagneticFieldVolume mVol, TVector3 pos);

    Bool_t FieldLoaded() { return GetNumberOfVolumes() == fMagneticFieldVolumes.size(); }

    MagneticFieldVolume* GetMagneticVolume(Int_t id) {
        if (!FieldLoaded()) LoadMagneticVolumes();
        if (fMagneticFieldVolumes.size() > id)
            return &fMagneticFieldVolumes[id];
        else {
            RESTError << "TRestAxionMagneticField::GetMagneticVolume. Id outside limits!" << RESTendl;
            return NULL;
        }
    }

   public:
    void LoadMagneticVolumes();

    Bool_t IsFieldConstant(Int_t id) {
        if (GetMagneticVolume(id)) return GetMagneticVolume(id)->field.size() == 0;
        return true;
    }

    Int_t GetNumberOfVolumes() { return fPositions.size(); }

    Bool_t CheckOverlaps();

    std::vector<TVector3> GetVolumeBoundaries(TVector3 pos, TVector3 dir, Int_t id = 0);
    std::vector<TVector3> GetFieldBoundaries(TVector3 pos, TVector3 dir, Double_t precision = 0,
                                             Int_t id = 0);

    TVector3 GetMagneticField(Double_t x, Double_t y, Double_t z);
    TVector3 GetMagneticField(TVector3 pos, Bool_t showWarning = true);

    Double_t GetPhotonMass(Double_t x, Double_t y, Double_t z, Double_t en);
    Double_t GetPhotonMass(TVector3 pos, Double_t en);
    Double_t GetPhotonMass(Int_t id, Double_t en);

    Double_t GetPhotonAbsorptionLength(Double_t x, Double_t y, Double_t z, Double_t en);
    Double_t GetPhotonAbsorptionLength(TVector3 pos, Double_t en);
    Double_t GetPhotonAbsorptionLength(Int_t id, Double_t en);

    Int_t GetVolumeIndex(TVector3 pos);
    Bool_t IsInside(TVector3 pos);

    TVector3 GetVolumePosition(Int_t id);
    TVector3 GetVolumeCenter(Int_t id);

    Double_t GetTransversalComponent(TVector3 position, TVector3 direction);

    std::vector<Double_t> GetTransversalComponentAlongPath(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);

    TVector3 GetFieldAverageTransverseVector(TVector3 from, TVector3 to, Double_t dl = 10., Int_t Nmax = 0);

    TCanvas* DrawHistogram(TString projection, TString Bcomp, Int_t volIndex = -1, Double_t step = -1,
                           TString style = "COLZ0", Double_t depth = -100010.0);

    TCanvas* DrawTracks(TVector3 vanishingPoint, Int_t divisions, Int_t volId = 0);

    void PrintMetadata();

    TRestAxionMagneticField();
    TRestAxionMagneticField(const char* cfgFileName, std::string name = "");
    ~TRestAxionMagneticField();

    ClassDef(TRestAxionMagneticField, 3);
};
#endif