TRestDetectorGas.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 RestCore_TRestDetectorGas
#define RestCore_TRestDetectorGas
 
#include <TApplication.h>
#include <TArrayI.h>
#include <TAxis.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TNamed.h>
#include <TROOT.h>
#include <TRestMetadata.h>
#include <TString.h>
#include <TSystem.h>
#include <TVector3.h>
 
#include <cstdlib>
#include <fstream>
#include <iostream>
 
#include "TRestDetectorDriftVolume.h"
 
#if defined USE_Garfield
#include <ComponentConstant.hh>
#include <GeometrySimple.hh>
#include <MediumMagboltz.hh>
#include <Sensor.hh>
#include <SolidBox.hh>
#include <TrackHeed.hh>
using namespace Garfield;
#else
class MediumMagboltz;
#endif
 
const int RESTGAS_ERROR = -1;
const int RESTGAS_INTITIALIZED = 0;
const int RESTGAS_CFG_LOADED = 1;
const int RESTGAS_GASFILE_LOADED = 2;
 
class TRestDetectorGas : public TRestDetectorDriftVolume {
   private:
    MediumMagboltz* fGasMedium;  
 
    Int_t fStatus;         // Used to define the status of the gas : RESTGAS_ERROR,
                           // RESTGAS_INTITIALIZED, RESTGAS_CFG_LOADED,
                           // RESTGAS_GASFILE_LOADED
    TString fGasFilename;  // The filename of the Magboltz gas file.
    Int_t fNofGases;       // Number of different elements composing the gas mixture
 
    Int_t fNCollisions;           // Number of collisions used in the Magboltz calculation.
    Double_t fMaxElectronEnergy;  // Maximum electron energy, in eV, used in
                                  // Magboltz gas calculation.
 
    std::vector<TString> fGasComponentName;       // A std::string std::vector storing the names
                                                  // of each of the gas components
    std::vector<Double_t> fGasComponentFraction;  // A double std::vector storing the fraction values of
                                                  // each of the gas components
 
    Int_t fEnodes;   // Number of electric field nodes used in the gas calculation.
    Double_t fEmax;  // Minimum value of the electric field used for the gas
                     // calculation.
    Double_t fEmin;  // Maximum value of the electric field used for the gas
                     // calculation.
 
    std::vector<Double_t> fEFields;  // The electric field nodes as calculated by
                                     // Garfield::MediumMagboltz.
    std::vector<Double_t> fBFields;  // The magnetic field nodes as calculated by
                                     // Garfield::MediumMagboltz.
    std::vector<Double_t> fAngles;   // The field angles as calculated by Garfield::MediumMagboltz.
 
    TString fGDMLMaterialRef;  // The corresponding material reference name in
                               // GDML description
 
    bool fGasGeneration;  
 
    TString fGasOutputPath;  
 
    TString fGasServer;  
 
    TString fGasFileContent;  // Used for saving the gasFile into a root file
 
    Bool_t fTest = false;  
 
    void InitFromConfigFile() override;
    std::string ConstructFilename();
 
    void AddGasComponent(const std::string& gasName, Double_t fraction);
 
    void GenerateGasFile();
 
    void UploadGasToServer(std::string gasFilename);
 
    Double_t GetDriftVelocity(Double_t E);
    Double_t GetLongitudinalDiffusion(Double_t E);
    Double_t GetTransversalDiffusion(Double_t E);
    Double_t GetTownsendCoefficient(Double_t E);
    Double_t GetAttachmentCoefficient(Double_t E);
 
   public:
    TRestDetectorGas();
    TRestDetectorGas(const char* configFilename, const std::string& name = "", bool gasGeneration = false,
                     bool test = false);
    ~TRestDetectorGas();
 
    inline void EnableGasGeneration() { fGasGeneration = true; }
 
    inline bool GasFileGenerationEnabled() const { return fGasGeneration; }
 
    inline bool GasFileLoaded() const { return fStatus == RESTGAS_GASFILE_LOADED; }
 
    void Initialize() override;
 
    void LoadGasFile();
 
    std::string FindGasFile(std::string name);
 
    void CalcGarField(double Emin, double Emax, int n);
 
    Int_t Write(const char* name = nullptr, Int_t option = 0, Int_t bufsize = 0) override;
 
    void InitFromRootFile() override;
 
    inline Double_t GetMaxElectronEnergy() const { return fMaxElectronEnergy; }
 
    inline Int_t GetNofGases() const { return fNofGases; }
 
    inline TString GetGasComponentName(Int_t n) {
        if (n >= GetNofGases()) {
            std::cout << "REST WARNING. Gas name component n=" << n << " requested. But only "
                      << GetNofGases() << " component(s) in the mixture." << std::endl;
            return "";
        }
        return fGasComponentName[n];
    }
 
    TString GetGasMixture();
 
    Double_t GetDriftVelocity() override {
        if (fElectricField == 0) {
            RESTWarning << "TRestDetectorGas::GetDriftVelocity. Warning fElectricField is zero!" << RESTendl;
            RESTWarning << " - Use: TRestDetectorGas::SetElectricField( field[V/mm] ) to set the field value"
                        << RESTendl;
        }
        return GetDriftVelocity(fElectricField * units("V/cm")) / units("cm/us");
    }  // in standard unit mm/us
 
    Double_t GetLongitudinalDiffusion() override {
        if (fElectricField == 0) {
            RESTWarning << "TRestDetectorGas::GetLongitudinalDiffusion. Warning fElectricField is zero!"
                        << RESTendl;
            RESTWarning << " - Use: TRestDetectorGas::SetElectricField( field[V/mm] ) to set the field value"
                        << RESTendl;
        }
        return GetLongitudinalDiffusion(fElectricField * units("V/cm"));
    }
 
    Double_t GetTransversalDiffusion() override {
        if (fElectricField == 0) {
            RESTWarning << "TRestDetectorGas::GetTransversalDiffusion. Warning fElectricField is zero!"
                        << RESTendl;
            RESTWarning << " - Use: TRestDetectorGas::SetElectricField( field[V/mm] ) to set the field value"
                        << RESTendl;
        }
        return GetTransversalDiffusion(fElectricField * units("V/cm"));
    }
 
    Double_t GetTownsendCoefficient() override {
        if (fElectricField == 0) {
            RESTWarning << "TRestDetectorGas::GetTownsendCoefficient. Warning fElectricField is zero!"
                        << RESTendl;
            RESTWarning << " - Use: TRestDetectorGas::SetElectricField( field[V/mm] ) to set the field value"
                        << RESTendl;
        }
        return GetTownsendCoefficient(fElectricField * units("V/cm"));
    }
 
    Double_t GetAttachmentCoefficient() override {
        if (fElectricField == 0) {
            RESTWarning << "TRestDetectorGas::GetAttachmentCoefficient. Warning fElectricField is zero!"
                        << RESTendl;
            RESTWarning << " - Use: TRestDetectorGas::SetElectricField( field[V/mm] ) to set the field value"
                        << RESTendl;
        }
        return GetAttachmentCoefficient(fElectricField * units("V/cm"));
    }
 
    void GetGasWorkFunctionOld();
 
    inline Double_t GetGasComponentFraction(Int_t n) {
        if (n >= GetNofGases()) {
            std::cout << "REST WARNING. Gas fraction for component n=" << n << " requested. But only "
                      << GetNofGases() << " component(s) in the mixture." << std::endl;
            return 0.;
        }
 
        return fGasComponentFraction[n];
    }
 
    Double_t GetGasFanoFactor() const;
 
    Double_t GetGasWorkFunction() const;
 
    inline MediumMagboltz* GetGasMedium() const { return fGasMedium; };
 
    inline TString GetGDMLMaterialRef() const { return fGDMLMaterialRef; };
 
    void SetPressure(Double_t pressure) override;
    void SetTemperature(Double_t temperature) override;
 
    inline void SetMaxElectronEnergy(Double_t energy) { fMaxElectronEnergy = energy; }
 
    void PlotDriftVelocity(Double_t eMin, Double_t eMax, Int_t nSteps);
    void PlotLongitudinalDiffusion(Double_t eMin, Double_t eMax, Int_t nSteps);
    void PlotTransversalDiffusion(Double_t eMin, Double_t eMax, Int_t nSteps);
    void PlotTownsendCoefficient(Double_t eMin, Double_t eMax, Int_t nSteps);
    void PrintGasInfo();
    inline void PrintGasFileContent() { std::cout << fGasFileContent << std::endl; };
 
    void PrintMetadata() override { PrintGasInfo(); }
 
    ClassDefOverride(TRestDetectorGas, 3);  // Gas Parameters
};
 
#endif