Graduation-Project-Gpp/mdt.cpp

#include <TApplication.h>
#include <TCanvas.h>
#include <TROOT.h>

#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/TrackSrim.hh"
#include "Garfield/ViewCell.hh"
#include "Garfield/ViewDrift.hh"
#include "Garfield/ViewSignal.hh"

#include <cstdlib>
#include <ctime>
#include <fstream>
#include <iostream>

#define f_PLOT 1
#define f_PRINT 2

using namespace Garfield;

const double rWire = 2e-3;
const double rTube = 2.54 / 2;
const double vWire = 1900.;
const double vTube = 0.;

const double tStep = 1;
const int nbins = 10000;
const int N = 1;

#if f_PRINT
std::ofstream outfile;
#endif

bool readTransferFunction(Sensor& sensor, std::string filename) {
    std::ifstream infile;
    infile.open(filename, std::ios::in);
    if (!infile) {
        std::cerr << "Could not read response function.\n";
        return false;
    }
    std::vector<double> times;
    std::vector<double> values;
    while (!infile.eof()) {
        double t = 0., f = 0.;
        infile >> t >> f;
        if (infile.eof() || infile.fail()) break;
        times.push_back(1.e3 * t);
        values.push_back(f);
    }
    infile.close();
    sensor.SetTransferFunction(times, values);
    return true;
}

double riseTime(Sensor& sensor) {
    int tl = -1, tr;
    double max = -1;
    double signal[nbins];

    for (int i = 0; i < nbins; i++) {
        signal[i] = -sensor.GetSignal("s", i);
        if (signal[i] > max) max = signal[i];
    }
    for (int i = 0; i < nbins; i++)
        if (signal[i] > max * 0.1) {
            tl = i;
            break;
        }
    for (int i = tl; i < nbins; i++)
        if (signal[i] > max * 0.6) {
            tr = i;
            break;
        }
    return 1.e-9 * (tr - tl) * tStep;
}

void randPos(double* x, double* y, double* z) {
    double a1 = 2 * M_PI * (rand() / (double)RAND_MAX);
    double a2 = (rand() / (double)RAND_MAX) * rTube;
    *x = cos(a1) * a2 + rWire;
    *y = sin(a1) * a2 + rWire;
    *z = (2 * (rand() / (double)RAND_MAX) - 1) * 10;
}

void randDir(double* x, double* y, double* z) {
    double a1 = 2 * M_PI * (rand() / (double)RAND_MAX);
    double a2 = acos(2 * (rand() / (double)RAND_MAX) - 1);
    *x = cos(a1) * cos(a2);
    *y = sin(a1) * cos(a2);
    *z = sin(a2);
}

std::tuple<int, double> driftElectron(TrackSrim track, DriftLineRKF drift) {
    int nc, ne = 0;
    double eDep = 0;
    double xc, yc, zc, tc, eDepc, extra;

    while (track.GetCluster(xc, yc, zc, tc, nc, eDepc, extra)) {
        ne += nc;
        eDep += eDepc;

        drift.SetIonSignalScalingFactor(nc);
        drift.DriftIon(xc, yc, zc, tc);
        drift.SetElectronSignalScalingFactor(nc);
        drift.DriftElectron(xc, yc, zc, tc);
#if f_PRINT
        outfile << xc << " " << yc << " " << zc << " " << tc << std::endl;
#endif
    }

    return {ne, eDep};
}

int main(int argc, char* argv[]) {
#if f_PRINT
    outfile.open("test.txt", std::ios::out);
#endif

    srand(time(NULL));
    TApplication app("app", &argc, argv);

    MediumMagboltz gas;
    gas.LoadGasFile("helium3_100atm_3000K.gas");
    gas.LoadIonMobility("IonMobility_He+_He.txt");

    ComponentAnalyticField cmp;
    cmp.SetMedium(&gas);

    cmp.AddWire(0., 0., 2 * rWire, vWire, "s", 30.);
    cmp.AddTube(rTube, vTube, 0, "t");
    cmp.AddReadout("s");

    Sensor sensor;
    sensor.AddComponent(&cmp);
    sensor.AddElectrode(&cmp, "s");
    sensor.SetTimeWindow(0, tStep, nbins);
    readTransferFunction(sensor, "mdt_CR_RC.txt");

    TrackSrim proton;
    proton.SetSensor(&sensor);
    proton.ReadFile("proton_in_he3(gas).txt");
    proton.SetKineticEnergy(0.573e6);

    TrackSrim tritium;
    tritium.SetSensor(&sensor);
    tritium.ReadFile("tritium_in_he3(gas).txt");
    tritium.SetKineticEnergy(0.191e6);

    TrackSrim alpha;
    alpha.SetSensor(&sensor);
    alpha.ReadFile("alpha_in_he3(gas).txt");

    DriftLineRKF drift;
    drift.SetSensor(&sensor);
    drift.SetMaximumStepSize(5.e-4);
    drift.EnableIonTail();
    drift.EnableSignalCalculation();

#if f_PLOT
    TCanvas* cD = new TCanvas("cD", "", 600, 600);
    ViewCell cellView;
    ViewDrift driftView;
    ViewSignal signalView;

    cellView.SetCanvas(cD);
    cellView.SetComponent(&cmp);
    driftView.SetCanvas(cD);
    drift.EnablePlotting(&driftView);
    proton.EnablePlotting(&driftView);
    tritium.EnablePlotting(&driftView);
    alpha.EnablePlotting(&driftView);
#endif

    double ek, rt, px, py, pz, dx, dy, dz;
    for (int i = 0; i < N; i++) {
#if f_PLOT
        driftView.Clear();
#endif
        sensor.ClearSignal();

        randPos(&px, &py, &pz);
        randDir(&dx, &dy, &dz);
        ek = 0.8e6 * (rand() / (double)RAND_MAX) + 0.2e6;

        // proton.SetKineticEnergy(ek);
        proton.NewTrack(px, py, pz, 0, dx, dy, dz);
        // tritium.SetKineticEnergy(ek);
        tritium.NewTrack(px, py, pz, 0, -dx, -dy, -dz);
        // alpha.SetKineticEnergy(ek);
        // alpha.NewTrack(px, py, pz, 0, dx, dy, dz);

        auto [ne, eDep] = driftElectron(proton, drift);
        auto [ne2, eDep2] = driftElectron(tritium, drift);

        // sensor.IntegrateSignals();
        sensor.ConvoluteSignals();
        rt = riseTime(sensor) * 1e9;

#if f_PLOT
        cD->Clear();
        cellView.Plot2d();
        driftView.Plot(true, false);
        cD->Update();

        signalView.SetSensor(&sensor);
        // signalView.PlotSignal("s", "ei");
        signalView.PlotSignal("s", "t");
#endif

#if f_PRINT > 1
        printf("Energy Deposit = %.3f, Rise Time = %.3f\n", eDep + eDep2, rt);
        printf("Position = (%.2f, %.2f), Direction = (%.2f, %.2f, %.2f)\n", px, py, dx, dy, dz);
#endif
#if f_PRINT
        // outfile << eDep + eDep2 << " " << rt << std::endl;
        outfile << sensor.GetSignal("s", 0);
        for (int i = 1; i < nbins; i++) outfile << " " << sensor.GetSignal("s", i);
        outfile << std::endl;
#endif
    }

    app.Run(kTRUE);
    return 0;
}