Graduation-Project-Gpp/proton.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>

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 = 20000;
const int N = 1000;

bool readTransferFunction(Sensor& sensor) {
    std::ifstream infile;
    infile.open("mdt_CR_RC.txt", std::ios::in);
    if (!infile) {
        std::cerr << "Could not read delta 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;
}

double rorate(double* x, double* y, double* z, double theta, int axis) {
    if (axis == 0) {
        double tmp = *y * cos(theta) - *z * sin(theta);
        *z = *y * sin(theta) + *z * cos(theta);
        *y = tmp;
    } else if (axis == 1) {
        double tmp = *x * cos(theta) - *z * sin(theta);
        *z = *x * sin(theta) + *z * cos(theta);
        *x = tmp;
    } else {
        double tmp = *x * cos(theta) - *y * sin(theta);
        *y = *x * sin(theta) + *y * cos(theta);
        *x = tmp;
    }
}

int main(int argc, char* argv[]) {
    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(-tStep / 2, tStep, nbins);
    readTransferFunction(sensor);

    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();

    // 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);

    int nc, ne, status;
    double eDep;
    double a1, a2, px, py, dx, dy, dz;
    double xc, yc, zc, tc, eDepc, extra;

    for (int i = 0; i < N; i++) {
        ne = 0;
        eDep = 0;
        // driftView.Clear();
        sensor.ClearSignal();

        a1 = 2 * M_PI * (rand() / (double)RAND_MAX);
        a2 = acos(2 * (rand() / (double)RAND_MAX) - 1);

        dx = cos(a1) * cos(a2);
        dy = sin(a1) * cos(a2);
        dz = sin(a2);

        a1 = 2 * M_PI * (rand() / (double)RAND_MAX);
        a2 = (rand() / (double)RAND_MAX) * rTube * 0.95;
        px = cos(a1) * a2 + rWire;
        py = sin(a1) * a2 + rWire;

        proton.SetKineticEnergy(0.6e6 * (rand() / (double)RAND_MAX) + 0.2e6);
        proton.NewTrack(px, py, 0, 0, dx, dy, dz);
        while (proton.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);
        }

        // tritium.NewTrack(px, py, 0, 0, -dx, -dy, -dz);
        // while (tritium.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);
        // }

        // alpha.NewTrack(rTrack, rTrack, 0, 0, -1, 1, 0);
        // while (alpha.GetCluster(xc, yc, zc, tc, nc, eDepc, extra)) {
        //     drift.SetIonSignalScalingFactor(nc);
        //     drift.DriftIon(xc, yc, zc, tc);
        //     drift.SetElectronSignalScalingFactor(nc);
        //     drift.DriftElectron(xc, yc, zc, tc);
        // }

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

        // // sensor.IntegrateSignals();
        // sensor.ConvoluteSignals();
        // signalView.SetSensor(&sensor);
        // // signalView.PlotSignal("s", "tei");
        // signalView.PlotSignal("s", "t");

        std::cout << eDep << " " << riseTime(sensor) * 1e9 << std::endl;
        // printf("Position = (%.2f, %.2f), ", px, py);
        // printf("Total Count = %d, Energy Deposit = %.3f, W = %.3f\n", ne, eDep, eDep / ne);
    }

    app.Run(kTRUE);
    return 0;
}