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 "TRestMesh.h"
#include "TVector3.h"
#include "TVectorD.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;  //<
 
    TVector3 fReMap = TVector3(0, 0, 0);  //<
 
    std::vector<MagneticFieldVolume> fMagneticFieldVolumes;  
 
    TVector3 fTrackStart = TVector3(0, 0, 0);  
 
    TVector3 fTrackDirection = TVector3(0, 0, 0);  
 
    Double_t fTrackLength = 0;  
 
    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 (id >= 0 && fMagneticFieldVolumes.size() > (unsigned int)id)
            return &fMagneticFieldVolumes[id];
        else {
            RESTError << "TRestAxionMagneticField::GetMagneticVolume. Id outside limits!" << RESTendl;
            return NULL;
        }
    }
 
   public:
    void LoadMagneticVolumes();
 
    void ReMap(const size_t& n, const TVector3& newMapSize);
 
    void SetTrack(const TVector3& position, const TVector3& direction);
 
    Double_t GetTrackLength() const { return fTrackLength; }
    TVector3 GetTrackStart() const { return fTrackStart; }
    TVector3 GetTrackDirection() const { return fTrackDirection; }
 
    Bool_t IsFieldConstant(Int_t id) {
        if (GetMagneticVolume(id)) return GetMagneticVolume(id)->field.size() == 0;
        return true;
    }
 
    unsigned int 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);
    Double_t GetTransversalComponentInParametricTrack(Double_t x);
 
    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