SciFiMatG4_1layer1mm/SciFiSim/src/Parameters.cc

// Written by Peter Stromberger
// based on work of Mirco Deckenhoff
// modified by Bastian Rössler




#include "Analysis.hh"
#include "Parameters.hh"

#include <fstream>
#include <vector>



Parameters* Parameters::singleton = 0;

Parameters::Parameters()
{

    //// Read parameter values and set corresponding variables ////

    G4cout << ">>> Constructor of Parameter class called <<<" << G4endl;

    //char parameterFileName[40] = "parameterFiles/parameters.dat";
    
#if defined __APPLE__
    char parameterFileName[200] = "/Users/basti/Programming/LHCb/SciFiSim-Xcode/SimulationData/parameterFiles/parameters.dat";
    #pragma message("Compiling for MacOSX")
#else
    char parameterFileName[40] = "parameterFiles/parameters.dat";
    #pragma message("Compiling for Linux")
#endif
    std::ifstream parameterFile;
    parameterFile.open(parameterFileName);

    if(parameterFile.good())
    {
        parameterFile >> randomSeed; // random seed for random engine
        parameterFile >> randomNumber; // random number for detector strip positioning
        randomNumber = randomNumber/32767.0; // max random number in bourne shell random number generator
        parameterFile >> fibreLength; // length of fibre in meter
        parameterFile >> semiAxisZ; // semi axis of fibre in z in millimeter
        parameterFile >> semiAxisY; // semi axis of fibre in y in millimeter
        parameterFile >> triggerX; // x-size of Trigger in mm
        parameterFile >> triggerY; // y-size of Trigger in mm
        parameterFile >> triggerZ; // z-size of Trigger in mm
        parameterFile >> triggerXPos; // x-position of Trigger in millimeter
        parameterFile >> triggerZPos; // z-position of Trigger in millimeter
        parameterFile >> probabilityOfPhotonLossAtSurface; // probability that a photon is killed when reaching fibre surface
      
        if(probabilityOfPhotonLossAtSurface<0 || probabilityOfPhotonLossAtSurface>1)
          probabilityOfPhotonLossAtSurface = 1;
      
        parameterFile >> placeMirror; // place a mirror at fibre end
        parameterFile >> mirrorReflectivity; // reflectivity of mirror at fibre end
        parameterFile >> detectorMaterial; // place a mirror at fibre end

        parameterFile.ignore(256,'\n');
        parameterFile.peek();
        parameterFile.getline(emissionSpectrumFileName,256);
        
        G4cout << "Emission FIle: " << emissionSpectrumFileName << G4endl;
        
        std::ifstream emissionFile;
        emissionFile.open(emissionSpectrumFileName);
      
        if(emissionFile.good())
        {
            std::vector<G4double> energyVector;
            std::vector<G4double> intensityVector;
            G4double readValue1;
            G4double readValue2;
            numberOfEnergies = 0;
          
            while(emissionFile >> readValue1)
            {
                if(emissionFile >> readValue2)
                {
                    energyVector.push_back(readValue1);
                    intensityVector.push_back(readValue2);
                    numberOfEnergies ++;
                }
            }
      
            Energy = new G4double[numberOfEnergies];
            Intensity = new G4double[numberOfEnergies];
          
            for(int i=0; i<numberOfEnergies; i++)
            {
                Energy[i] = energyVector[i];
                Intensity[i] = intensityVector[i];
            }
        }
        else
        {
            G4cout << "Could not read emission spectrum file: " << emissionSpectrumFileName <<" !" << G4endl;
        }
      
        emissionFile.close();
        parameterFile >> numberOfInterpolatedPoints; // number of points to be interpolated within each emission spectrum interval

        parameterFile.ignore(256,'\n');
        parameterFile.peek();
        parameterFile.getline(wlsAbsSpectrumFileName,256);
        std::ifstream wlsAbsFile;
        wlsAbsFile.open(wlsAbsSpectrumFileName);
    
        if(wlsAbsFile.good())
        {
            std::vector<G4double> wlsAbsEnergyVector;
            std::vector<G4double> wlsAbsLengthVector;
            G4double readValue1;
            G4double readValue2;
            numberOfWlsAbsEnergies = 0;
          
            while(wlsAbsFile >> readValue1)
            {
                if(wlsAbsFile >> readValue2)
                {
                    wlsAbsEnergyVector.push_back(readValue1);
                    wlsAbsLengthVector.push_back(readValue2);
                    numberOfWlsAbsEnergies ++;
                }
            }
      
            WlsAbsEnergy = new G4double[numberOfWlsAbsEnergies];
            WlsAbsLength = new G4double[numberOfWlsAbsEnergies];
          
            for(int i=0; i<numberOfWlsAbsEnergies; i++)
            {
                WlsAbsEnergy[i] = wlsAbsEnergyVector[i];
                WlsAbsLength[i] = wlsAbsLengthVector[i];
            }
        }
        else
        {
            G4cout << "Could not read wls absorption spectrum file: " << wlsAbsSpectrumFileName <<" !" << G4endl;
        }
      
        wlsAbsFile.close();

        parameterFile.getline(wlsEmissionSpectrumFileName,256);
        std::ifstream wlsEmissionFile;
        wlsEmissionFile.open(wlsEmissionSpectrumFileName);
    
        if(wlsEmissionFile.good())
        {
            std::vector<G4double> wlsEmissionEnergyVector;
            std::vector<G4double> wlsEmissionIntensityVector;
            G4double readValue1;
            G4double readValue2;
            numberOfWlsEmissionEnergies = 0;
            
            while(wlsEmissionFile >> readValue1)
            {
                if(wlsEmissionFile >> readValue2)
                {
                    wlsEmissionEnergyVector.push_back(readValue1);
                    wlsEmissionIntensityVector.push_back(readValue2);
                    numberOfWlsEmissionEnergies ++;
                }
            }
      
            WlsEmissionEnergy = new G4double[numberOfWlsEmissionEnergies];
            WlsEmissionIntensity = new G4double[numberOfWlsEmissionEnergies];
      
            for(int i=0; i<numberOfWlsEmissionEnergies; i++)
            {
                WlsEmissionEnergy[i] = wlsEmissionEnergyVector[i];
                WlsEmissionIntensity[i] = wlsEmissionIntensityVector[i];
            }
        }
        else
        {
            G4cout << "Could not read wls emission spectrum file: " << wlsEmissionSpectrumFileName <<" !" << G4endl;
        }
        
        wlsEmissionFile.close();

        parameterFile >> scintillationYield; //scintillation yield in photons per keV
        parameterFile >> resolutionScale; // width of gaussian to generate photon number
        parameterFile >> decayTimeFast; // fast decay time of excited states in ns
        parameterFile >> decayTimeSlow; // slow decay time of excited states in ns
        parameterFile >> yieldRatio; // ratio of fast component and total scintillation yield
        parameterFile >> birksConstant; // Birk's constant in mm/MeV

        parameterFile >> wlsDecayTime; // decay time of wls excited states in ns

        parameterFile.ignore(256,'\n');
        parameterFile.peek();
        parameterFile.getline(refractiveIndexVacuum,256);
        parameterFile.getline(refractiveIndexCore,256);
        parameterFile.getline(refractiveIndexClad1,256);
        parameterFile.getline(refractiveIndexClad2,256);

        parameterFile.getline(absorptionCore,512);
        parameterFile.getline(absorptionClad1,512);
        parameterFile.getline(absorptionClad2,512);

        parameterFile.getline(absorptionFromIrradiationCore,512);
        parameterFile.getline(absorptionFromIrradiationClad1,512);
        parameterFile.getline(absorptionFromIrradiationClad2,512);

        parameterFile.getline(sectionsFileName,80);
        std::ifstream sectionsFile;
        sectionsFile.open(sectionsFileName);
    
        if(sectionsFile.good())
            sectionsFile >> numberOfSections;
        else
            G4cout << "Could not read sections file: " << sectionsFileName <<" !" << G4endl;
    
        sectionsFile.close();

        parameterFile.getline(rayleighCore,512);
        parameterFile.getline(rayleighClad1,512);
        parameterFile.getline(rayleighClad2,512);
    }
    else
    {
        G4cout << "Could not read parameter file: " << parameterFileName <<" !" << G4endl;
    }

    parameterFile.close();



  //// Open file to store parameters ////

    sprintf(parameterOutputFileName,"%s.parameters",Analysis::GetInstance()->FileName());

    std::ofstream parameterOutputFile;
    parameterOutputFile.open(parameterOutputFileName);

    parameterOutputFile << "Length of fibre: ";
    parameterOutputFile << fibreLength << " m \n";
    parameterOutputFile << "Semi axis of fibre in z: ";
    parameterOutputFile << semiAxisZ << " mm \n";
    parameterOutputFile << "Semi axis of fibre in y: ";
    parameterOutputFile << semiAxisY << " mm \n";
    parameterOutputFile << "Probability to lose photons at fibre surface: ";
    parameterOutputFile << probabilityOfPhotonLossAtSurface << "\n";
    parameterOutputFile << "Mirror placement at fibre end: ";
    parameterOutputFile << placeMirror << "\n";
    parameterOutputFile << "Reflectivity of mirror at fibre end: ";
    parameterOutputFile << mirrorReflectivity << "\n";
    parameterOutputFile << "Detector material vacuum/polystyrene (0/1): ";
    parameterOutputFile << detectorMaterial << "\n";
    parameterOutputFile << "Used emission spectrum: \"";
    parameterOutputFile << emissionSpectrumFileName << "\"\n";
    parameterOutputFile << "Number of energies: ";
    parameterOutputFile << numberOfEnergies << "\n";
    parameterOutputFile << "Energy / eV \t Intensity\n";
  
    for(int i=0; i<numberOfEnergies; i++)
        parameterOutputFile << Energy[i] << "\t" << Intensity[i] <<"\n";
  
    parameterOutputFile << "Number of interpolated points per emission spectrum interval: ";
    parameterOutputFile << numberOfInterpolatedPoints << "\n";
    parameterOutputFile << "Used wls absorption spectrum: \"";
    parameterOutputFile << wlsAbsSpectrumFileName << "\"\n";
    parameterOutputFile << "Number of wls absorption energies: ";
    parameterOutputFile << numberOfWlsAbsEnergies << "\n";
    parameterOutputFile << "Energy / eV \t WLS Absorption Length / m\n";
  
    for(int i=0; i<numberOfWlsAbsEnergies; i++)
        parameterOutputFile << WlsAbsEnergy[i] << "\t" << WlsAbsLength[i] <<"\n";
  
    parameterOutputFile << "Used wls emission spectrum: \"";
    parameterOutputFile << wlsEmissionSpectrumFileName << "\"\n";
    parameterOutputFile << "Number of wls emission energies: ";
    parameterOutputFile << numberOfWlsEmissionEnergies << "\n";
    parameterOutputFile << "Energy / eV \t WLS emission intensity\n";
    
    for(int i=0; i<numberOfWlsEmissionEnergies; i++)
        parameterOutputFile << WlsEmissionEnergy[i] << "\t" << WlsEmissionIntensity[i] <<"\n";
  
    parameterOutputFile << "Formulae to calculate refractive indices in\n";
    parameterOutputFile << " - vacuum: \"" << refractiveIndexVacuum << "\"\n";
    parameterOutputFile << " - core: \"" << refractiveIndexCore << "\"\n";
    parameterOutputFile << " - inner cladding: \"" << refractiveIndexClad1 << "\"\n";
    parameterOutputFile << " - outer cladding: \"" << refractiveIndexClad2 << "\"\n";
    parameterOutputFile << " with x in nm.\n";
    parameterOutputFile << "Scintillation yield: " << scintillationYield << "/keV\n";
    parameterOutputFile << "Resolution scale: " << resolutionScale << "\n";
    parameterOutputFile << "Fast decay time: " << decayTimeFast << " ns\n";
    parameterOutputFile << "Slow decay time: " << decayTimeSlow << " ns\n";
    parameterOutputFile << "Yield ratio: " << yieldRatio << "\n";
    parameterOutputFile << "Birks constant: " << birksConstant << " mm/MeV\n";
    parameterOutputFile << "WLS decay time: " << wlsDecayTime << " ns\n";
    parameterOutputFile << "Formulae to calculate absorption [1/m] in\n";
    parameterOutputFile << " - core: \"" << absorptionCore << "\"\n";
    parameterOutputFile << " - inner cladding: \"" << absorptionClad1 << "\"\n";
    parameterOutputFile << " - outer cladding: \"" << absorptionClad2 << "\"\n";
    parameterOutputFile << " with x in nm.\n";
    parameterOutputFile << "Formulae to calculate absorption due to irradiation [1/m] in\n";
    parameterOutputFile << " - core: \"" << absorptionFromIrradiationCore << "\"\n";
    parameterOutputFile << " - inner cladding: \"" << absorptionFromIrradiationClad1 << "\"\n";
    parameterOutputFile << " - outer cladding: \"" << absorptionFromIrradiationClad2 << "\"\n";
    parameterOutputFile << " with x in nm and y in kGy.\n";
    parameterOutputFile << "Used sections file: \"";
    parameterOutputFile << sectionsFileName << "\"\n";
    parameterOutputFile << "Number of sections: ";
    parameterOutputFile << numberOfSections << "\n";
    parameterOutputFile << "Formulae to calculate Rayleigh scattering [1/m] in\n";
    parameterOutputFile << " - core: " << rayleighCore << "\n";
    parameterOutputFile << " - inner cladding: " << rayleighClad1 << "\n";
    parameterOutputFile << " - outer cladding: " << rayleighClad2 << "\n";

    parameterOutputFile.close();

    G4cout << "Simulation parameters written to \"" << parameterOutputFileName << "\"" << G4endl;
}




Parameters::~Parameters() {
    if(numberOfEnergies != 0)
    {
        delete[] Energy;
        delete[] Intensity;
    }
    if(numberOfWlsAbsEnergies != 0)
    {
        delete[] WlsAbsEnergy;
        delete[] WlsAbsLength;
    }
    if(numberOfWlsEmissionEnergies != 0)
    {
        delete[] WlsEmissionEnergy;
        delete[] WlsEmissionIntensity;
    }
}





char* Parameters::ParameterOutputFileName()
{
    return parameterOutputFileName;
}

G4double Parameters::RandomSeed()
{
    return randomSeed;
}

G4double Parameters::RandomNumber()
{
    return randomNumber;
}

G4double Parameters::FibreLength()
{
    return fibreLength;
}


G4double Parameters::SemiAxisZ()
{
    return semiAxisZ;
}


G4double Parameters::SemiAxisY()
{
    return semiAxisY;
}


G4double Parameters::TriggerX()
{
    return triggerX;
}

G4double Parameters::TriggerY()
{
	return triggerY;
}

G4double Parameters::TriggerZ()
{
	return triggerZ;
}

G4double Parameters::TriggerXPos()
{
	return triggerXPos;
}

G4double Parameters::TriggerZPos()
{
	return triggerZPos;
}

G4double Parameters::ProbabilityOfPhotonLossAtSurface()
{
    return probabilityOfPhotonLossAtSurface;
}


G4bool Parameters::PlaceMirror()
{
    return placeMirror;
}


G4double Parameters::MirrorReflectivity()
{
    return mirrorReflectivity;
}


G4bool Parameters::DetectorMaterial()
{
    return detectorMaterial;
}


char* Parameters::EmissionSpectrumFileName()
{
    return emissionSpectrumFileName;
}


G4int Parameters::NumberOfEnergies()
{
    return numberOfEnergies;
}


G4int Parameters::NumberOfInterpolatedPoints()
{
    return numberOfInterpolatedPoints;
}


char* Parameters::WlsAbsSpectrumFileName()
{
    return wlsAbsSpectrumFileName;
}


G4int Parameters::NumberOfWlsAbsEnergies()
{
    return numberOfWlsAbsEnergies;
}


char* Parameters::WlsEmissionSpectrumFileName()
{
    return wlsEmissionSpectrumFileName;
}


G4int Parameters::NumberOfWlsEmissionEnergies()
{
    return numberOfWlsEmissionEnergies;
}


G4double Parameters::ScintillationYield()
{
    return scintillationYield;
}


G4double Parameters::ResolutionScale()
{
    return resolutionScale;
}
 

G4double Parameters::DecayTimeFast()
{
  return decayTimeFast;
}
 

G4double Parameters::DecayTimeSlow()
{
    return decayTimeSlow;
}
 

G4double Parameters::YieldRatio()
{
    return yieldRatio;
}
 

G4double Parameters::BirksConstant()
{
    return birksConstant;
}


G4double Parameters::WlsDecayTime()
{
    return wlsDecayTime;
}
 

char* Parameters::RefractiveIndexVacuum()
{
  return refractiveIndexVacuum;
}


char* Parameters::RefractiveIndexCore()
{
    return refractiveIndexCore;
}


char* Parameters::RefractiveIndexClad1()
{
    return refractiveIndexClad1;
}


char* Parameters::RefractiveIndexClad2()
{
    return refractiveIndexClad2;
}


char* Parameters::AbsorptionCore()
{
    return absorptionCore;
}


char* Parameters::AbsorptionClad1()
{
    return absorptionClad1;
}


char* Parameters::AbsorptionClad2()
{
    return absorptionClad2;
}

 
char* Parameters::AbsorptionFromIrradiationCore()
{
    return absorptionFromIrradiationCore;
}


char* Parameters::AbsorptionFromIrradiationClad1()
{
    return absorptionFromIrradiationClad1;
}


char* Parameters::AbsorptionFromIrradiationClad2()
{
    return absorptionFromIrradiationClad1;
}


char* Parameters::SectionsFileName()
{
    return sectionsFileName;
}


G4int Parameters::NumberOfSections()
{
    return numberOfSections;
}


char* Parameters::RayleighCore()
{
    return rayleighCore;
}


char* Parameters::RayleighClad1()
{
    return rayleighClad1;
}


char* Parameters::RayleighClad2()
{
    return rayleighClad2;
}


const G4double  Parameters::hcPERe = 1.239842e-6;  // unit: V*m