InductionGar2.cxx

Go to the documentation of this file.

#include "CgemDigitizerSvc/InductionGar2.h"
 
#include "CLHEP/Units/PhysicalConstants.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/DataSvc.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "G4DigiManager.hh"
#include "G4ios.hh"
#include "Randomize.hh"
 
#include "CLHEP/Units/PhysicalConstants.h"
 
#include "TRandom.h"
#include "TTree.h"
#include "TSystem.h"
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
 
//#define INDUCTIONGAR_DEBUG_RECTANGLE2_FAST_TEST
//uncomment the previous line to check correctness of results of rectangle2_fast.
//it will be slow.
 
//#define INDUCTIONGAR_DEBUG_CONVOLUTION_TBRANCH_FFT_TEST
//uncomment the previous line to check correctness of results of Convolution_Tbranch_fft and Convolution_Ebranch_fft
//it will be very slow.
 
//#define INDUCTIONGAR_DEBUG_CONV1PERGRID_FFT_TEST
//uncomment the previous line to check equality of 2 conv1PerGrid.
 
const static int electrons_select_method_threhold = 0; //always use FFT conv here.
 
const static int _nbins = 2000;
typedef std::vector<int> Vint;
using namespace std;
const static double zeros[_nbins] = { 0};
static void compareMap(const std::map<int, double> mapQ1[2][2], const std::map<int, double> mapQ2[2][2]);
static void DrawToFile(const double *h, const char *filename);
static void checkMemorySize();
 
 
template <typename T>
static bool compareVector(const vector<T> &v1, const vector<T> &v2, const char *description, double threshold = 0, bool noOutput = false) {
    if (threshold < 0) threshold = 0;
    //static const double threshold=1e-5;
    //# if std::numeric_limits<T>::is_signed
    bool found_mismatch      = false;
    const vector<T> *vLarger = &v1;
    char vLargerName[]       = "v1";
    if (noOutput){
        if (v1.size() != v2.size()) return false;
        for (size_t i = 0; i < std::min(v1.size(), v2.size()); i++) {
            if (max(v1[i], v2[i]) - min(v1[i], v2[i]) > threshold) return false;
        }
    }else{
        if (v1.size() != v2.size()) {
            {
                if (not found_mismatch) {
                    cout << "compareVector:" << description << ":found mismatch:" << '\n';
                    found_mismatch = true;
                }
            }
            cout << "\t size not match, v1:" << v1.size() << " v2:" << v2.size() << "\n";
            if (v1.size() < v2.size()) {
                vLarger        = &v2;
                vLargerName[1] = '2';
            }
        }
 
        for (size_t i = 0; i < std::min(v1.size(), v2.size()); i++) {
            if (max(v1[i], v2[i]) - min(v1[i], v2[i]) > threshold) {
                {
                    if (not found_mismatch) {
                        cout << "compareVector:" << description << ":found mismatch:" << '\n';
                        found_mismatch = true;
                    }
                }
                cout << "\t at i=" << i << " not match, v1[i] = " << v1[i] << " v2[i] = " << v2[i] << '\n';
            }
        }
 
        for (size_t i = min(v1.size(), v2.size()); i < max(v1.size(), v2.size()); i++) {
            {
                if (not found_mismatch) {
                    cout << "compareVector:" << description << ":found mismatch:" << '\n';
                    found_mismatch = true;
                }
            }
            cout << "\t at i=" << i << ", only has " << vLargerName << "[i] = " << (*vLarger)[i] << '\n';
        }
    }
    return !found_mismatch;
}
 
template <typename T>
static bool compareArray(const T *v1, const T *v2, size_t length, const char *description, double threshold = 0, bool noOutput=false) {
    if (threshold < 0) threshold = 0;
    //static const double threshold=1e-5;
    //# if std::numeric_limits<T>::is_signed
    bool found_mismatch = false;
 
    for (size_t i = 0; i < length; i++) {
        if (max(v1[i], v2[i]) - min(v1[i], v2[i]) > threshold) {
            {
                if (noOutput) return false;
                if (not found_mismatch) {
                    cout << "compareArray:" << description << ":found mismatch:" << '\n';
                    found_mismatch = true;
                }
            }
            cout << "\t at i=" << i << " not match, v1[i] = " << v1[i] << " v2[i] = " << v2[i] << '\n';
        }
    }
    return !found_mismatch;
}
 
//-----------------------Convolution of response function---------------------------------
 
static double rectangle2(double t, double *x, double *y, int Npx) {
 
    double I = y[0];
    int id   = -1;
    for (int i = 0; i < Npx - 1; i++) {
        if (t >= x[i] && t < x[i + 1]) {
            id = i;
            break;
        }
    }
    if (id != -1) {
        I = y[id] + (y[id + 1] - y[id]) * (t - x[id]) / (x[id + 1] - x[id]);
    }
 
    return I;
}
 
/**
 * @brief calc estimate using linear interpolation of curve {x,y} , at point t;  mimic rectangle2() but faster version
 *      assuming x: from 0 to 1000 and 2000 points.
 * @param t 
 * @param x 
 * @param y 
 * @param Npx 
 * @return double 
 */
static double rectangle2_fast(double t, double *x, const double *y, int Npx) {
 
    const int nbinsx      = 2000;
    const double xmin     = 0;
    const double xmax     = 1000;
    const double binWidth = (xmax - xmin) / nbinsx;
    double I              = y[0];
    int id                = -1;
 
    if (t > xmin + binWidth / 2 && t < xmax - binWidth / 2) {
        id = (int)floor((t - binWidth / 2 - xmin) / binWidth);
    }
#ifdef INDUCTIONGAR_DEBUG_RECTANGLE2_FAST_TEST // checking results of new and old version. for speed it should not run
    int id2 = -1;
    for (int i = 0; i < Npx - 1; i++) {
        if (t >= x[i] && t < x[i + 1]) {
            id2 = i;
            break;
        } // find the closist id, to make t=x[id]; find root.
    }
    if (id2 != id) {
        cout << "CgemDigitizerSvc::InductionGar2::rectangle2_fast: 2 method calced id doesnot correspond, "
             << "id new : " << id2 << " id old : " << id << " t : " << t << " x at id old " << x[id2] << '\n';
    }
#endif
    if (id != -1) {
        I = y[id] + (y[id + 1] - y[id]) * (t - x[id]) / (x[id + 1] - x[id]);
    } // the tangent of y(i),x(i) at i=id, I is on tangent and at x=t;
 
    return I;
}
 
static double rectangle2(double t, std::vector<double> &x, const double *y, int Npx) {
 
    double I = y[0];
    int id   = -1;
    for (int i = 0; i < Npx - 1; i++) {
        if (t >= x[i] && t < x[i + 1]) {
            id = i;
            break;
        }
    }
    if (id != -1) {
        I = y[id] + (y[id + 1] - y[id]) * (t - x[id]) / (x[id + 1] - x[id]);
    }
 
    return I;
}
 
static double T_branch2(double t) {
    double result = 0.;
    t             = t - 10.;
    if (t > 0) {
        result = 2000. * (0.00181928 * exp(-t / 3.) - 0.0147059 * exp(-t / 20.) + 0.0128866 * exp(-t / 100.));
    }
    return result;
}
 
static double E_branch2(double t) {
    double result = 0.;
    t             = t - 5;
    if (t > 0) {
        result = 0.000627357 * (1358.7 * exp(-t * 0.0385647) * t + 1358.7 * exp(-t * 0.0114353) * t + 100164. * exp(-t * 0.0385647) - 100164. * exp(-0.0114353 * t));
    }
    return result;
}
 
static TH1F Convolution_Tbranch(const double *Input_Curr_plot_001, double xmin, double xmax, int Npx = 2000) // return TIGER output
{
 
    TH1F h_signalT("h_signalT", "", Npx, xmin, xmax);
    std::vector<double> Input_Time_plot_001(Npx);
    for (int i = 0; i < Npx; i++) Input_Time_plot_001[i] = h_signalT.GetBinCenter(i + 1);
 
    double xmin_f1(0.), xmax_f1(1000);
    int nStep_f1(200);
    double step_f1 = (xmax_f1 - xmin_f1) / nStep_f1;
    double x_f1(0.0);
 
    std::vector<double> x(Npx);
    std::vector<double> y_001(Npx);
 
    for (int i = 0; i < Npx; i++) {
        x[i]           = xmin + (xmax - xmin) / Npx * (i + 0.5);
        double fun_001 = 0.0;
 
        for (int j = 0; j < nStep_f1; j++) {
            x_f1 = xmin_f1 + step_f1 * (j + 0.5);
            fun_001 += rectangle2(x_f1, Input_Time_plot_001, Input_Curr_plot_001, Npx) * T_branch2(x[i] - x_f1) * step_f1;
        }
        y_001[i] = fun_001;
    }
 
    for (int init = 0; init < Npx; init++) {
        h_signalT.SetBinContent(init + 1, -y_001[init]);
    }
 
    return h_signalT;
}
 
static TH1F Convolution_Tbranch(const double *Input_Curr_plot_001) // return TIGER output
{
    double xmin(0), xmax(1000);
    const int Npx(2000);
 
    TH1F h_signalT("h_signalT", "", Npx, xmin, xmax);
    double Input_Time_plot_001[Npx];
    for (int i = 0; i < Npx; i++) Input_Time_plot_001[i] = h_signalT.GetBinCenter(i + 1);
 
    double xmin_f1(0.), xmax_f1(1000);
    int nStep_f1(200);
    double step_f1 = (xmax_f1 - xmin_f1) / nStep_f1;
    double x_f1(0.0);
 
    double x[Npx];
    double y_001[Npx];
 
    for (int i = 0; i < Npx; i++) {
        x[i]           = xmin + (xmax - xmin) / Npx * (i + 0.5);
        double fun_001 = 0.0;
 
        for (int j = 0; j < nStep_f1; j++) {
            x_f1 = xmin_f1 + step_f1 * (j + 0.5);
            fun_001 += rectangle2_fast(x_f1, Input_Time_plot_001, Input_Curr_plot_001, Npx) * T_branch2(x[i] - x_f1) * step_f1;
        }
        y_001[i] = fun_001;
    }
 
    for (int init = 0; init < Npx; init++) {
        h_signalT.SetBinContent(init + 1, -y_001[init]);
    }
 
    return h_signalT;
}
 
static TH1F Convolution_Ebranch(const double *Input_Curr_plot_001, double xmin, double xmax, int Npx = 2000) // return TIGER output
{
    TH1F h_signalE("h_signalE", "", Npx, xmin, xmax);
    std::vector<double> Input_Time_plot_001(Npx);
    for (int i = 0; i < Npx; i++) Input_Time_plot_001[i] = h_signalE.GetBinCenter(i + 1);
 
    double xmin_f1(0.), xmax_f1(1000);
    int nStep_f1(200);
    double step_f1 = (xmax_f1 - xmin_f1) / nStep_f1;
    double x_f1(0.0);
 
    std::vector<double> x(Npx);
    std::vector<double> y_001(Npx);
 
    for (int i = 0; i < Npx; i++) {
        x[i]           = xmin + (xmax - xmin) / Npx * (i + 0.5);
        double fun_001 = 0.0;
 
        for (int j = 0; j < nStep_f1; j++) {
            x_f1 = xmin_f1 + step_f1 * (j + 0.5);
            fun_001 += rectangle2(x_f1, Input_Time_plot_001, Input_Curr_plot_001, Npx) * E_branch2(x[i] - x_f1) * step_f1;
        }
        y_001[i] = fun_001;
    }
 
    for (int init = 0; init < Npx; init++) {
        h_signalE.SetBinContent(init + 1, -y_001[init]);
    }
 
    return h_signalE;
}
 
static TH1F Convolution_Ebranch(const double *Input_Curr_plot_001) // return TIGER output
{
    double xmin(0), xmax(1000);
    const int Npx(2000);
 
    TH1F h_signalE("h_signalE", "", Npx, xmin, xmax);
    double Input_Time_plot_001[Npx];
    for (int i = 0; i < Npx; i++) Input_Time_plot_001[i] = h_signalE.GetBinCenter(i + 1);
 
    double xmin_f1(0.), xmax_f1(1000);
    int nStep_f1(200);
    double step_f1 = (xmax_f1 - xmin_f1) / nStep_f1;
    double x_f1(0.0);
 
    double x[Npx];
    double y_001[Npx];
 
    for (int i = 0; i < Npx; i++) {
        x[i]           = xmin + (xmax - xmin) / Npx * (i + 0.5);
        double fun_001 = 0.0;
 
        for (int j = 0; j < nStep_f1; j++) {
            x_f1 = xmin_f1 + step_f1 * (j + 0.5);
            fun_001 += rectangle2_fast(x_f1, Input_Time_plot_001, Input_Curr_plot_001, Npx) * E_branch2(x[i] - x_f1) * step_f1;
        }
        y_001[i] = fun_001;
    }
 
    for (int init = 0; init < Npx; init++) {
        h_signalE.SetBinContent(init + 1, -y_001[init]);
    }
 
    return h_signalE;
}
 
#ifdef INDUCTIONGAR_DEBUG_CONVOLUTION_TBRANCH_FFT_TEST
/*
 * @brief works exactly like Convolution_Tbranch, but works on every point rather than 1 in 10. 
 * 
 * @param Input_Curr_plot_001 
 * @return TH1F 
 */
static TH1F Convolution_Tbranch_noskip(const double *Input_Curr_plot_001) // return TIGER output
{
    double xmin(0), xmax(1000);
    const int Npx(2000);
 
    TH1F h_signalT("h_signalT", "", Npx, xmin, xmax);
    double Input_Time_plot_001[Npx];
    for (int i = 0; i < Npx; i++) Input_Time_plot_001[i] = h_signalT.GetBinCenter(i + 1);
 
    double xmin_f1(0.), xmax_f1(1000);
    int nStep_f1(2000);
    double step_f1 = (xmax_f1 - xmin_f1) / nStep_f1;
    double x_f1(0.0);
 
    double x[Npx];
    double y_001[Npx];
 
    for (int i = 0; i < Npx; i++) {
        x[i]           = xmin + (xmax - xmin) / Npx * (i + 0.5);
        double fun_001 = 0.0;
 
        for (int j = 0; j < nStep_f1; j++) {
            x_f1 = xmin_f1 + step_f1 * (j + 0.5);
            fun_001 += rectangle2_fast(x_f1, Input_Time_plot_001, Input_Curr_plot_001, Npx) * T_branch2(x[i] - x_f1) * step_f1;
        }
        y_001[i] = fun_001;
    }
 
    for (int init = 0; init < Npx; init++) {
        h_signalT.SetBinContent(init + 1, -y_001[init]);
    }
 
    return h_signalT;
}
 
static TH1F Convolution_Ebranch_noskip(const double *Input_Curr_plot_001) // return TIGER output
{
    double xmin(0), xmax(1000);
    const int Npx(2000);
 
    TH1F h_signalE("h_signalE", "", Npx, xmin, xmax);
    double Input_Time_plot_001[Npx];
    for (int i = 0; i < Npx; i++) Input_Time_plot_001[i] = h_signalE.GetBinCenter(i + 1);
 
    double xmin_f1(0.), xmax_f1(1000);
    int nStep_f1(2000);
    double step_f1 = (xmax_f1 - xmin_f1) / nStep_f1;
    double x_f1(0.0);
 
    double x[Npx];
    double y_001[Npx];
 
    for (int i = 0; i < Npx; i++) {
        x[i]           = xmin + (xmax - xmin) / Npx * (i + 0.5);
        double fun_001 = 0.0;
 
        for (int j = 0; j < nStep_f1; j++) {
            x_f1 = xmin_f1 + step_f1 * (j + 0.5);
            fun_001 += rectangle2_fast(x_f1, Input_Time_plot_001, Input_Curr_plot_001, Npx) * E_branch2(x[i] - x_f1) * step_f1;
        }
        y_001[i] = fun_001;
    }
 
    for (int init = 0; init < Npx; init++) {
        h_signalE.SetBinContent(init + 1, -y_001[init]);
    }
 
    return h_signalE;
}
#endif
CTF2::CTF2() {
    double xmin(0), xmax(1000);
    const int Npx(2000);
    TBranch.resize(Npx); // TBranch at x<0 is 0
    double *p_TBranch = &(TBranch.front());
    for (int i = 0; i < Npx; i++) {
        double x   = xmin + (xmax - xmin) * i / Npx;
        TBranch[i] = T_branch2(x);
    } // set initial TBranch values.
    // basically, a N and a (2N-1) needs (3N-2) bins to avoid overlap.
    conv = new ConvolveWithConst(p_TBranch, Npx, Npx);
    //TBranch_fft.resize(Npx*3-2);
 
    //conv.convolve
}
CTF2::~CTF2() {
    delete conv;
}
/**
 * @brief return convolve of Input_Curr_plot_001 and T_Branch2, a function; tries to mimic Convolution_Tbranch
 * has a option to test correctness.
 * 
 * @param Input_Curr_plot_001 
 * @return TH1D 
 */
TH1D CTF2::Convolution_Tbranch_fft(const double *Input_Curr_plot_001) // return TIGER output; note it returns DOUBLE typed hist, rather than original float.
{
    double xmin(0), xmax(1000);
    const int Npx(2000);
 
    TH1D h_signal("h_signal", "", Npx, xmin, xmax);
    double *p_signal = h_signal.GetArray() + 1;                        //considering underflow.
    this->conv->convolve(p_signal, Input_Curr_plot_001, 0, Npx, -0.5); // the '-' before 0.5: from original function.
 
#ifdef INDUCTIONGAR_DEBUG_CONVOLUTION_TBRANCH_FFT_TEST
    const double dx_abs_threhold = 1e-6; // think about the accurancy of float; warn if  |x| > threhold
    const double dx_ref_threhold = 1e-6; //think about the accurancy of float; warn if |dx/x| > threhold
    TH1F h_signal_ref            = Convolution_Tbranch_noskip(Input_Curr_plot_001);
    bool found_mismatch          = false;
    for (int i = 1; i <= Npx; i++) {
        if (abs((h_signal[i] - h_signal_ref[i]) / h_signal_ref[i]) > dx_ref_threhold and abs(h_signal[i] - h_signal_ref[i]) > dx_abs_threhold) {
            if (not found_mismatch) {
                std::cout << "CgemDigitizerSvc::InductionGar2::Convolution_Tbranch_fft found results not match:  ";
                found_mismatch = true;
            }
            std::cout << "ref: " << h_signal_ref[i] << "; this: " << h_signal[i] << "; at i=" << i << ";threhold abs=" << dx_abs_threhold << "; ref=" << dx_ref_threhold << '\n';
        }
    }
#endif
 
    return h_signal;
}
 
CEF2::CEF2() {
    double xmin(0), xmax(1000);
    const int Npx(2000);
    EBranch.resize(Npx); // EBranch at x<0 is 0
    double *p_EBranch = &(EBranch.front());
    for (int i = 0; i < Npx; i++) {
        double x   = xmin + (xmax - xmin) * i / Npx;
        EBranch[i] = E_branch2(x);
    } // set initial EBranch values.
    // basically, a N and a (2N-1) needs (3N-2) bins to avoid overlap.
    conv = new ConvolveWithConst(p_EBranch, Npx, Npx);
    //EBranch_fft.resize(Npx*3-2);
 
    //conv.convolve
}
CEF2::~CEF2() {
    delete conv;
}
/**
 * @brief return convolve of Input_Curr_plot_001 and T_Branch2, a function; tries to mimic Convolution_Ebranch
 * has a option to test correctness.
 * 
 * @param Input_Curr_plot_001 
 * @return TH1D 
 */
TH1D CEF2::Convolution_Ebranch_fft(const double *Input_Curr_plot_001) // return TIGER output; note it returns DOUBLE typed hist, rather than original float.
{
    double xmin(0), xmax(1000);
    const int Npx(2000);
 
    TH1D h_signal("h_signalE", "", Npx, xmin, xmax);
    double *p_signal = h_signal.GetArray() + 1; //considering underflow.
    this->conv->convolve(p_signal, Input_Curr_plot_001, 0, Npx, -0.5);
 
#ifdef INDUCTIONGAR_DEBUG_CONVOLUTION_TBRANCH_FFT_TEST
    const double dx_abs_threhold = 1e-6; // think about the accurancy of float; warn if  |x| > threhold
    const double dx_ref_threhold = 1e-6; //think about the accurancy of float; warn if |dx/x| > threhold
    TH1F h_signal_ref            = Convolution_Ebranch_noskip(Input_Curr_plot_001);
    bool found_mismatch          = false;
    for (int i = 1; i <= Npx; i++) {
        if (abs((h_signal[i] - h_signal_ref[i]) / h_signal_ref[i]) > dx_ref_threhold and abs(h_signal[i] - h_signal_ref[i]) > dx_abs_threhold) {
            if (not found_mismatch) {
                std::cout << "CgemDigitizerSvc::InductionGar2::Convolution_Ebranch_fft found results not match:  ";
                found_mismatch = true;
            }
            std::cout << "ref: " << h_signal_ref[i] << "; this: " << h_signal[i] << "; at i=" << i << ";threhold abs=" << dx_abs_threhold << "; ref=" << dx_ref_threhold << '\n';
        }
    }
#endif
 
    return h_signal;
}
 
InductionGar2::InductionGar2() {
    string testPath = getenv("CGEMDIGITIZERSVCROOT");
    m_LUTFilePath   = "/bes3fs/cgemCosmic/data/CGEM_cosmic_look_up_table_from_10_to_17.root";
    m_V2EfineFile   = testPath + "testFIle/V2Efine.root";
    sample_delay    = 162.5;
    storeFlag       = false;
    m_saturation    = true;
 
 
    m_QinGausSigma[0]=2;
    m_QinGausSigma[1]=2;
}
 
InductionGar2::~InductionGar2() {
}
 
void InductionGar2::init(ICgemGeomSvc *geomSvc, double magConfig) {
    m_geomSvc   = geomSvc;
    m_magConfig = magConfig;
    std::cout << "InductionGar2 sampling time : " << sample_delay << std::endl;
 
    string filePath_ratio = getenv("CGEMDIGITIZERSVCROOT");
    string fileName;
    if (0 == m_magConfig)
        fileName = filePath_ratio + "/dat/par_InductionGar_0T.txt";
    else
        fileName = filePath_ratio + "/dat/par_InductionGar_1T.txt";
    ifstream fin(fileName.c_str(), ios::in);
    cout << "InductionGar2 fileName: " << fileName << endl;
 
    string strline;
    string strpar;
    if (!fin.is_open()) {
        cout << "ERROR: can not open file " << fileName << endl;
    }
 
    for (int m = 0; m < 3; m++) {
        for (int l = 0; l < 5; l++) {
            for (int k = 0; k < 5; k++) {
                for (int i = 0; i < 4; i++) {
                    fin >> RatioX[m][l][k][i];
                }
                for (int j = 0; j < 3; j++) {
                    fin >> RatioV[m][l][k][j];
                }
            }
        }
    }
    fin.close();
 
    string fileName1 = filePath_ratio + "/dat/VQ_relation.txt";
    ifstream VQin(fileName1.c_str(), ios::in);
    VQin >> VQ_slope >> VQ_const;
    VQin.close();
    //cout<<VQ_slope<<"    "<<VQ_const<<endl;
 
    string filePath = getenv("CGEMDIGITIZERSVCROOT");
    for (int i = 0; i < 3; i++) {         // layer
        for (int j = 0; j < 5; j++) {     // Grid-x
            for (int k = 0; k < 5; k++) { // Grid-v
                char temp1[1];
                sprintf(temp1, "%i", (i + 1));
                char temp2[2];
                sprintf(temp2, "%i", ((j + 1) * 10 + k + 1));
                if (0 == m_magConfig)
                    fileName = filePath + "/dat/InducedCurrent_Root_0T/layer" + temp1 + "_" + temp2 + ".root";
                else
                    fileName = filePath + "/dat/InducedCurrent_Root_1T/layer" + temp1 + "_" + temp2 + ".root";
                TFile *file_00   = TFile::Open(fileName.c_str(), "read");
                TH1 *readThis_x3 = 0;
                TH1 *readThis_x4 = 0;
                TH1 *readThis_x5 = 0;
                TH1 *readThis_x6 = 0;
                TH1 *readThis_v5 = 0;
                TH1 *readThis_v6 = 0;
                TH1 *readThis_v7 = 0;
                file_00->GetObject("readThis_x3", readThis_x3);
                file_00->GetObject("readThis_x4", readThis_x4);
                file_00->GetObject("readThis_v5", readThis_v5);
                file_00->GetObject("readThis_x5", readThis_x5);
                file_00->GetObject("readThis_v6", readThis_v6);
                file_00->GetObject("readThis_x6", readThis_x6);
                file_00->GetObject("readThis_v7", readThis_v7);
 
                double scaleX=m_ScaleSignalX;
                double scaleV=(1-scaleX*(RatioX[i][j][k][0]+RatioX[i][j][k][1]+RatioX[i][j][k][2]+RatioX[i][j][k][3]))/(RatioV[i][j][k][0]+RatioV[i][j][k][1]+RatioV[i][j][k][2]);
                for (int init = 0; init < 100; init++) {
                    SignalX[i][j][k][0][init] = scaleX*readThis_x3->GetBinContent(init + 1);
                    SignalX[i][j][k][1][init] = scaleX*readThis_x4->GetBinContent(init + 1);
                    SignalX[i][j][k][2][init] = scaleX*readThis_x5->GetBinContent(init + 1);
                    SignalX[i][j][k][3][init] = scaleX*readThis_x6->GetBinContent(init + 1);
                    SignalV[i][j][k][0][init] = scaleV*readThis_v5->GetBinContent(init + 1);
                    SignalV[i][j][k][1][init] = scaleV*readThis_v6->GetBinContent(init + 1);
                    SignalV[i][j][k][2][init] = scaleV*readThis_v7->GetBinContent(init + 1);
                }
                double temp[7][200];
                for (int init = 0; init < 100; init++) { // since it was used like 2 bins eventually.
                    temp[0][init * 2]     = SignalX[i][j][k][0][init];
                    temp[1][init * 2]     = SignalX[i][j][k][1][init];
                    temp[2][init * 2]     = SignalX[i][j][k][2][init];
                    temp[3][init * 2]     = SignalX[i][j][k][3][init];
                    temp[4][init * 2]     = SignalV[i][j][k][0][init];
                    temp[5][init * 2]     = SignalV[i][j][k][1][init];
                    temp[6][init * 2]     = SignalV[i][j][k][2][init];
                    temp[0][init * 2 + 1] = SignalX[i][j][k][0][init];
                    temp[1][init * 2 + 1] = SignalX[i][j][k][1][init];
                    temp[2][init * 2 + 1] = SignalX[i][j][k][2][init];
                    temp[3][init * 2 + 1] = SignalX[i][j][k][3][init];
                    temp[4][init * 2 + 1] = SignalV[i][j][k][0][init];
                    temp[5][init * 2 + 1] = SignalV[i][j][k][1][init];
                    temp[6][init * 2 + 1] = SignalV[i][j][k][2][init];
                }
                // since we actually used only index 1~80.
                Signal_ConvolverX[i][j][k][0].init(temp[0] + 2, 160, 2000);
                Signal_ConvolverX[i][j][k][1].init(temp[1] + 2, 160, 2000);
                Signal_ConvolverX[i][j][k][2].init(temp[2] + 2, 160, 2000);
                Signal_ConvolverX[i][j][k][3].init(temp[3] + 2, 160, 2000);
                Signal_ConvolverV[i][j][k][0].init(temp[4] + 2, 160, 2000);
                Signal_ConvolverV[i][j][k][1].init(temp[5] + 2, 160, 2000);
                Signal_ConvolverV[i][j][k][2].init(temp[6] + 2, 160, 2000);
            }
        }
    }
 
    TFile *LUTFile = TFile::Open(m_LUTFilePath.c_str(), "read");
    TTree *tree    = (TTree *)LUTFile->Get("tree");
    Float_t V_th_T, V_th_E, QDC_a, QDC_b, QDC_saturation;
    Int_t layer, sheet, strip_v, strip_x;
    tree->SetBranchAddress("strip_x_boss", &strip_x);
    tree->SetBranchAddress("strip_v_boss", &strip_v);
    tree->SetBranchAddress("layer", &layer);
    tree->SetBranchAddress("sheet", &sheet);
    tree->SetBranchAddress("calib_QDC_slope", &QDC_a);
    tree->SetBranchAddress("calib_QDC_const", &QDC_b);
    tree->SetBranchAddress("calib_QDC_saturation", &QDC_saturation);
    tree->SetBranchAddress("v_thr_T_mV", &V_th_T);
    tree->SetBranchAddress("v_thr_E_mV", &V_th_E);
 
    double thresholdMin_X_T = 999;
    double thresholdMin_V_T = 999;
    double thresholdMin_X_Q = 999;
    double thresholdMin_V_Q = 999;
    for (int i = 0; i < tree->GetEntries(); i++) {
        tree->GetEntry(i);
        if (strip_x != -1) {
            T_thr_V[layer][sheet][0][strip_x] = V_th_T;
            if (V_th_T > 0 && V_th_T < thresholdMin_X_T) thresholdMin_X_T = V_th_T;
            //if(V_th_T<0) T_thr_V[layer][sheet][0][strip_x] =12.;
            E_thr_V[layer][sheet][0][strip_x] = V_th_E;
            if (V_th_E > 0 && V_th_E < thresholdMin_X_Q) thresholdMin_X_Q = V_th_E;
            //if(V_th_E<0) E_thr_V[layer][sheet][0][strip_x] =12.;
            QDC_slope[layer][sheet][0][strip_x]   = QDC_a;
            QDC_const[layer][sheet][0][strip_x]   = QDC_b;
            Qsaturation[layer][sheet][0][strip_x] = QDC_saturation;
        }
        if (strip_v != -1) {
            T_thr_V[layer][sheet][1][strip_v] = V_th_T;
            if (V_th_T > 0 && V_th_T < thresholdMin_V_T) thresholdMin_V_T = V_th_T;
            //if(V_th_T<0) T_thr_V[layer][sheet][1][strip_v] =12.;
            E_thr_V[layer][sheet][1][strip_v] = V_th_E;
            if (V_th_E > 0 && V_th_E < thresholdMin_V_Q) thresholdMin_V_Q = V_th_E;
            //if(V_th_E<0) E_thr_V[layer][sheet][1][strip_v] =12.;
            QDC_slope[layer][sheet][1][strip_v]   = QDC_a;
            QDC_const[layer][sheet][1][strip_v]   = QDC_b;
            Qsaturation[layer][sheet][1][strip_v] = QDC_saturation;
        }
    }
    for (int i = 0; i < tree->GetEntries(); i++) {
        tree->GetEntry(i);
        if (strip_x != -1) {
            if (V_th_T <= 0) T_thr_V[layer][sheet][0][strip_x] = thresholdMin_X_T;
            if (V_th_E <= 0) E_thr_V[layer][sheet][0][strip_x] = thresholdMin_X_Q;
        }
        if (strip_v != -1) {
            if (V_th_T <= 0) T_thr_V[layer][sheet][1][strip_v] = thresholdMin_V_T;
            if (V_th_E <= 0) E_thr_V[layer][sheet][1][strip_v] = thresholdMin_V_Q;
        }
    }
 
    //3rd layer
    for (int i = 0; i < 832; i++) {
        T_thr_V[2][0][0][i]     = 12.;
        E_thr_V[2][0][0][i]     = 12.;
        QDC_slope[2][0][0][i]   = -10.47;
        QDC_const[2][0][0][i]   = 482.3;
        Qsaturation[2][0][0][i] = 45.;
        T_thr_V[2][1][0][i]     = 10.;
        E_thr_V[2][1][0][i]     = 10.;
        QDC_slope[2][1][0][i]   = -10.47;
        QDC_const[2][1][0][i]   = 482.3;
        Qsaturation[2][1][0][i] = 45.;
    }
 
    for (int i = 0; i < 1395; i++) {
        T_thr_V[2][0][1][i]     = 12.;
        E_thr_V[2][0][1][i]     = 12.;
        QDC_slope[2][0][1][i]   = -10.47;
        QDC_const[2][0][1][i]   = 482.3;
        Qsaturation[2][0][1][i] = 45.;
        T_thr_V[2][1][1][i]     = 10.;
        E_thr_V[2][1][1][i]     = 10.;
        QDC_slope[2][1][1][i]   = -10.47;
        QDC_const[2][1][1][i]   = 482.3;
        Qsaturation[2][1][1][i] = 45.;
    }
 
    LUTFile->Close();
    /*
       TFile *V_Efine_File = new TFile(m_V2EfineFile.c_str());
       TTree *tree_V2E = (TTree*)V_Efine_File->Get("tree");
       Double_t V_ref,V2E_slope,V2E_const;
       Int_t layer_,sheet_,stripv,stripx;
       tree_V2E->SetBranchAddress("strip_x", &stripx);
       tree_V2E->SetBranchAddress("strip_v", &stripv);
       tree_V2E->SetBranchAddress("layer", &layer_);
       tree_V2E->SetBranchAddress("sheet", &sheet_);
       tree_V2E->SetBranchAddress("V2E_slope", &V2E_slope);
       tree_V2E->SetBranchAddress("V2E_const", &V2E_const);
       tree_V2E->SetBranchAddress("V_ref", &V_ref);
 
 
       for(int i=0;i<tree_V2E->GetEntries();i++){
       tree_V2E->GetEntry(i);
       if(stripx!=-1){
       Vref[layer_][sheet_][0][stripx] = V_ref;
       toE_slope[layer_][sheet_][0][stripx] = V2E_slope;
       toE_const[layer_][sheet_][0][stripx] = V2E_const;
       }
       if(stripv!=-1){
       Vref[layer_][sheet_][1][stripv] = V_ref;
       toE_slope[layer_][sheet_][1][stripv] = V2E_slope;
       toE_const[layer_][sheet_][1][stripv] = V2E_const;
       }
       }
 
       V_Efine_File->Close();
       */
 
    cout<<"InductionGar: "<<endl;
    cout<<"QinGausSigma="<<m_QinGausSigma[0]
        <<", "<<m_QinGausSigma[1]
        <<endl;
    //m_deadStrip[3][2]
    m_deadStrip[0][0][0].insert(40);//x_0_down
    m_deadStrip[0][0][0].insert(189);
    m_deadStrip[0][0][0].insert(322);
    m_deadStrip[0][0][0].insert(350);
    m_deadStrip[0][0][0].insert(457);//x_0_up
    m_deadStrip[0][0][1].insert(282);//v_0_down
    m_deadStrip[0][0][1].insert(467);
    m_deadStrip[0][0][1].insert(715);
    m_deadStrip[0][0][1].insert(773);
    m_deadStrip[0][0][1].insert(795);
    m_deadStrip[0][0][1].insert(803);
    m_deadStrip[1][0][1].insert(305);//v_1_down
    m_deadStrip[1][0][1].insert(307);
    m_deadStrip[1][0][1].insert(381);
    m_deadStrip[1][0][1].insert(443);
    m_deadStrip[1][0][1].insert(486);
    m_deadStrip[1][0][1].insert(550);
    m_deadStrip[1][0][1].insert(620);
    m_deadStrip[1][0][1].insert(631);
    m_deadStrip[1][0][1].insert(660);
    m_deadStrip[1][0][1].insert(681);
    m_deadStrip[1][1][1].insert(425);//v_1_up
    m_deadStrip[1][1][1].insert(455);
    m_deadStrip[1][1][1].insert(459);
    m_deadStrip[1][1][1].insert(461);
    m_deadStrip[1][1][1].insert(535);
    m_deadStrip[1][1][1].insert(536);
    m_deadStrip[1][1][1].insert(614);
    m_deadStrip[1][1][1].insert(672);
}
 
void InductionGar2::setMultiElectrons(int layer, const HitHistMap &hitmap, const PVectorXYZT *debug_ref_electrons) {
    bool bShowRef=false;
    if (debug_ref_electrons!=0)bShowRef=true;
 
    #pragma region
    m_xstripSheet.clear();
    m_xstripID.clear();
    m_vstripSheet.clear();
    m_vstripID.clear();
    m_xstripQ.clear();
    m_vstripQ.clear();
    m_xstripT_Branch.clear();
    m_vstripT_Branch.clear();
    m_xstripQ_Branch.clear();
    m_vstripQ_Branch.clear();
    m_xTfirst.clear();
    m_vTfirst.clear();
 
    //int m000_nXstrips;
    //int m000_nVstrips;
    std::vector<int> m000_xstripSheet;
    std::vector<int> m000_xstripID;
    std::vector<int> m000_vstripSheet;
    std::vector<int> m000_vstripID;
    std::vector<double> m000_xstripQ;
    std::vector<double> m000_vstripQ;
    std::vector<double> m000_xstripT_Branch;
    std::vector<double> m000_vstripT_Branch;
    std::vector<double> m000_xstripQ_Branch;
    std::vector<double> m000_vstripQ_Branch;
 
    //m000_xstripSheet.clear();
    //m000_xstripID.clear();
    //m000_vstripSheet.clear();
    //m000_vstripID.clear();
    //m000_xstripQ.clear();
    //m000_vstripQ.clear();
    //m000_xstripT_Branch.clear();
    //m000_vstripT_Branch.clear();
    //m000_xstripQ_Branch.clear();
    //m000_vstripQ_Branch.clear();
 
    CgemGeoLayer *CgemLayer;
    CgemGeoReadoutPlane *CgemReadoutPlane;
    CgemLayer  = m_geomSvc->getCgemLayer(layer);
    int NSheet = CgemLayer->getNumberOfSheet();
 
    Vint Save_Sheetx, Save_Sheetv;
    Vint Save_FinalstripIDx, Save_FinalstripIDv;
    Vint Save_Gridx, Save_Gridv;
    Vint Save_IDx, Save_IDv;
    Vint Save_Tx, Save_Tv; //CAUTION! precision: 1, or 2bins; bias -0.5
    //Save_Tx.clear();
    //Save_Tv.clear();
    //Save_IDx.clear();
    //Save_IDv.clear();
    //Save_Gridx.clear();
    //Save_Gridv.clear();
    //Save_Sheetx.clear();
    //Save_Sheetv.clear();
    //Save_FinalstripIDx.clear();
    //Save_FinalstripIDv.clear();
    for (int i = 0; i < 2; i++) {
        for (int k = 0; k < 2; k++) {
            m_mapQ[i][k].clear();
            //m_mapT[i][k].clear();
        }
    }
    #pragma endregion
    // this loop works on m_mapQ, Save_{Grid,T,ID}{xv} (represents info {grid, time, strip} of electron . looks all in hitmap), Save_Sheet{XV}( sheet_id*10000 + stripID )
    //std::cout << "nElectrons = " << nElectrons << std::endl;
    //std::cout << "size of HitHistMap: "<<hitmap.size() <<" x "<<hitmap.begin()->second.capacity() <<" x "<<sizeof(hitmap.begin()->second[0])<<endl;
 
    if (bShowRef) { // check if HitHistMap and electron list has same hit number.
        int Nelec                      = 0;
        HitHistMap::const_iterator itm = hitmap.begin();
        for (; itm != hitmap.end(); itm++) {
            const Vint &hist = itm->second;
            for (int num = 0; num < hist.size(); num++) {
                Nelec += hist[num];
            }
        }
        if (Nelec != debug_ref_electrons->x->size()) {
            cout << "oops! given different electrons! map gives " << Nelec << " while list gives" << debug_ref_electrons->x->size() << endl;
        }
    }
 
    
    if (bShowRef) {
        const vector<float> &x = *debug_ref_electrons->x;
        const vector<float> &y = *debug_ref_electrons->y;
        const vector<float> &z = *debug_ref_electrons->z;
        const vector<float> &t = *debug_ref_electrons->t;
        const int nElectrons   = x.size();
 
        for (int i = 0; i < nElectrons; i++) {
            //cout<<i<<endl;
            double phi_electron = atan2(y[i], x[i]);
            double z_electron   = z[i];
            G4ThreeVector pos(x[i], y[i], z[i]);
            for (int j = 0; j < NSheet; j++) {
                CgemReadoutPlane = m_geomSvc->getReadoutPlane(layer, j);
                if (CgemReadoutPlane->OnThePlane(phi_electron, z_electron)) {
                    int xID;
                    double distX = CgemReadoutPlane->getDist2ClosestXStripCenter(phi_electron, xID);
                    int vID;
                    double distV = CgemReadoutPlane->getDist2ClosestVStripCenter(pos, vID);
 
                    //std::cout << "-------------" << std::endl;
                    //std::cout << "xID, distX = " << xID << "," << distX << std::endl;
                    //std::cout << "vID, distV = " << vID << "," << distV << std::endl;
 
                    int grid_x = 1;
                    int grid_v = 1;
 
                    //
                    //--- Calculation of grid index-------------
                    //
                    double L_XPitch = CgemReadoutPlane->getXPitch();
                    double L_VPitch = CgemReadoutPlane->getVPitch();
                    grid_v          = floor(distX / L_XPitch * 5. + 2.5);
                    grid_x          = floor(distV / L_VPitch * 5. + 2.5);
                    //std::cout << i << " " << j << " " << grid_x << "  " << grid_v << std::endl;
                    //std::cout << "L_XPitch,L_VPitch = " << L_XPitch       << "," << L_VPitch       << std::endl;
                    //std::cout << "grid_x, grid_v = "    << grid_x         << "," << grid_v         << std::endl;
 
                    if (xID >= 0 && vID >= 0) {
                        map<int, double>::iterator itx = m_mapQ[j][0].find(xID);
                        map<int, double>::iterator itv = m_mapQ[j][1].find(vID);
                        if (RatioX[layer][grid_x][grid_v][2] != 0) {
                            if (itx == m_mapQ[j][0].end()) {
                                m_mapQ[j][0][xID] = RatioX[layer][grid_x][grid_v][2];
                                Save_Gridx.push_back(grid_x * 100 + grid_v * 10 + 2);
                                Save_IDx.push_back(j * 10000 + xID);
                                Save_Tx.push_back(t[i]);
                                //std::cout << i << " " << j << "x-strip : " << RatioX[layer][grid_x][grid_v][2] << std::endl;
                            } else {
                                double Q          = (itx->second) + RatioX[layer][grid_x][grid_v][2];
                                m_mapQ[j][0][xID] = Q;
                                Save_Gridx.push_back(grid_x * 100 + grid_v * 10 + 2);
                                Save_IDx.push_back(j * 10000 + xID);
                                Save_Tx.push_back(t[i]);
                                //std::cout << i << " " << j << "x-strip : Add " << Q << std::endl;
                            }
                            vector<int>::iterator result_Sheetx = find(Save_Sheetx.begin(), Save_Sheetx.end(), j * 10000 + xID);
                            if (result_Sheetx == Save_Sheetx.end()) {
                                Save_Sheetx.push_back(j * 10000 + xID);
                            }
                        }
                        if (RatioV[layer][grid_x][grid_v][1] != 0) {
                            if (itv == m_mapQ[j][1].end()) {
                                m_mapQ[j][1][vID] = RatioV[layer][grid_x][grid_v][1];
                                Save_Gridv.push_back(grid_x * 100 + grid_v * 10 + 1);
                                Save_IDv.push_back(j * 10000 + vID);
                                Save_Tv.push_back(t[i]);
                            } else {
                                double Q          = (itv->second) + RatioV[layer][grid_x][grid_v][1];
                                m_mapQ[j][1][vID] = Q;
                                Save_Gridv.push_back(grid_x * 100 + grid_v * 10 + 1);
                                Save_IDv.push_back(j * 10000 + vID);
                                Save_Tv.push_back(t[i]);
                            }
                            vector<int>::iterator result_Sheetv = find(Save_Sheetv.begin(), Save_Sheetv.end(), j * 10000 + vID);
                            if (result_Sheetv == Save_Sheetv.end()) {
                                Save_Sheetv.push_back(j * 10000 + vID);
                            }
                        }
                    }
 
                    if ((xID >= 0 && vID >= 0) && (xID + 1) < CgemReadoutPlane->getNXstrips()) {
                        map<int, double>::iterator itx = m_mapQ[j][0].find(xID + 1);
                        if (RatioX[layer][grid_x][grid_v][3] != 0) {
                            if (itx == m_mapQ[j][0].end()) {
                                m_mapQ[j][0][xID + 1] = RatioX[layer][grid_x][grid_v][3];
                                Save_Gridx.push_back(grid_x * 100 + grid_v * 10 + 3);
                                Save_IDx.push_back(j * 10000 + xID + 1);
                                Save_Tx.push_back(t[i]);
                                //std::cout << i << " " << j << "x-strip : " << RatioX[layer][grid_x][grid_v][3] << std::endl;
                            } else {
                                double Q              = (itx->second) + RatioX[layer][grid_x][grid_v][3];
                                m_mapQ[j][0][xID + 1] = Q;
                                Save_Gridx.push_back(grid_x * 100 + grid_v * 10 + 3);
                                Save_IDx.push_back(j * 10000 + xID + 1);
                                Save_Tx.push_back(t[i]);
                            }
 
                            vector<int>::iterator result_Sheetx = find(Save_Sheetx.begin(), Save_Sheetx.end(), j * 10000 + xID + 1);
                            if (result_Sheetx == Save_Sheetx.end()) {
                                Save_Sheetx.push_back(j * 10000 + xID + 1);
                            }
                        }
                    }
 
                    if ((xID >= 0 && vID >= 0) && (vID + 1) < CgemReadoutPlane->getNVstrips()) {
                        map<int, double>::iterator itv = m_mapQ[j][1].find(vID + 1);
                        if (RatioV[layer][grid_x][grid_v][2] != 0) {
                            if (itv == m_mapQ[j][1].end()) {
                                m_mapQ[j][1][vID + 1] = RatioV[layer][grid_x][grid_v][2];
                                Save_Gridv.push_back(grid_x * 100 + grid_v * 10 + 2);
                                Save_IDv.push_back(j * 10000 + vID + 1);
                                Save_Tv.push_back(t[i]);
                            } else {
                                double Q              = (itv->second) + RatioV[layer][grid_x][grid_v][2];
                                m_mapQ[j][1][vID + 1] = Q;
                                Save_Gridv.push_back(grid_x * 100 + grid_v * 10 + 2);
                                Save_IDv.push_back(j * 10000 + vID + 1);
                                Save_Tv.push_back(t[i]);
                            }
 
                            vector<int>::iterator result_Sheetv = find(Save_Sheetv.begin(), Save_Sheetv.end(), j * 10000 + vID + 1);
                            if (result_Sheetv == Save_Sheetv.end()) {
                                Save_Sheetv.push_back(j * 10000 + vID + 1);
                            }
                        }
                    }
 
                    if ((xID >= 0 && vID >= 0) && (xID - 2) >= 0) {
                        map<int, double>::iterator itx = m_mapQ[j][0].find(xID - 2);
                        if (RatioX[layer][grid_x][grid_v][0] != 0) {
                            if (itx == m_mapQ[j][0].end()) {
                                m_mapQ[j][0][xID - 2] = RatioX[layer][grid_x][grid_v][0];
                                Save_Gridx.push_back(grid_x * 100 + grid_v * 10 + 0);
                                Save_IDx.push_back(j * 10000 + xID - 2);
                                Save_Tx.push_back(t[i]);
                                //std::cout << i << " " << j << "x-strip : " << RatioX[layer][grid_x][grid_v][0] << std::endl;
                            } else {
                                double Q              = (itx->second) + RatioX[layer][grid_x][grid_v][0];
                                m_mapQ[j][0][xID - 2] = Q;
                                Save_Gridx.push_back(grid_x * 100 + grid_v * 10 + 0);
                                Save_IDx.push_back(j * 10000 + xID - 2);
                                Save_Tx.push_back(t[i]);
                            }
 
                            vector<int>::iterator result_Sheetx = find(Save_Sheetx.begin(), Save_Sheetx.end(), j * 10000 + xID - 2);
                            if (result_Sheetx == Save_Sheetx.end()) {
                                Save_Sheetx.push_back(j * 10000 + xID - 2);
                            }
                        }
                    }
 
                    if ((xID >= 0 && vID >= 0) && (xID - 1) >= 0) {
                        map<int, double>::iterator itx = m_mapQ[j][0].find(xID - 1);
                        if (RatioX[layer][grid_x][grid_v][1] != 0) {
                            if (itx == m_mapQ[j][0].end()) {
                                m_mapQ[j][0][xID - 1] = RatioX[layer][grid_x][grid_v][1];
                                Save_Gridx.push_back(grid_x * 100 + grid_v * 10 + 1);
                                Save_IDx.push_back(j * 10000 + xID - 1);
                                Save_Tx.push_back(t[i]);
                            } else {
                                double Q              = (itx->second) + RatioX[layer][grid_x][grid_v][1];
                                m_mapQ[j][0][xID - 1] = Q;
                                Save_Gridx.push_back(grid_x * 100 + grid_v * 10 + 1);
                                Save_IDx.push_back(j * 10000 + xID - 1);
                                Save_Tx.push_back(t[i]);
                            }
 
                            vector<int>::iterator result_Sheetx = find(Save_Sheetx.begin(), Save_Sheetx.end(), j * 10000 + xID - 1);
                            if (result_Sheetx == Save_Sheetx.end()) {
                                Save_Sheetx.push_back(j * 10000 + xID - 1);
                            }
                        }
                    }
 
                    if ((xID >= 0 && vID >= 0) && (vID - 1) >= 0) {
                        map<int, double>::iterator itv = m_mapQ[j][1].find(vID - 1);
                        if (RatioV[layer][grid_x][grid_v][0] != 0) {
                            if (itv == m_mapQ[j][1].end()) {
                                m_mapQ[j][1][vID - 1] = RatioV[layer][grid_x][grid_v][0];
                                Save_Gridv.push_back(grid_x * 100 + grid_v * 10 + 0);
                                Save_IDv.push_back(j * 10000 + vID - 1);
                                Save_Tv.push_back(t[i]);
                            } else {
                                double Q              = (itv->second) + RatioV[layer][grid_x][grid_v][0];
                                m_mapQ[j][1][vID - 1] = Q;
                                Save_Gridv.push_back(grid_x * 100 + grid_v * 10 + 0);
                                Save_IDv.push_back(j * 10000 + vID - 1);
                                Save_Tv.push_back(t[i]);
                            }
 
                            vector<int>::iterator result_Sheetv = find(Save_Sheetv.begin(), Save_Sheetv.end(), j * 10000 + vID - 1);
                            if (result_Sheetv == Save_Sheetv.end()) {
                                Save_Sheetv.push_back(j * 10000 + vID - 1);
                            }
                        }
                    }
                }
            }
        }
    }
 
    
 
    //std::cout << "loop electron finished " << std::endl;
 
    // the above loop comes to
    std::map<IndexGar,std::vector<double> > hitFFTMap; // size is 2* HalfPlus1 aka 2* (m_length/2+1) 
    
    if (bShowRef) {
        std::map<int, double> mapQ[2][2];
 
        calcMapQFromHitHistMap(hitmap, layer, mapQ);
        Vint Save_Sheetx2, Save_Sheetv2;
        calcSaveSheetXVFromMapQ(Save_Sheetx2, Save_Sheetv2, mapQ);
 
        compareMap(m_mapQ, mapQ);
 
        // TODO: check Save_Sheet* and Save_Sheet*2
        compareVector(Save_Sheetx, Save_Sheetx2, "Save_Sheetx", 1e-5);
        compareVector(Save_Sheetv, Save_Sheetv2, "Save_Sheetv", 1e-5);
    } else {
 
        calcMapQFromHitHistMap(hitmap, layer, m_mapQ);
        calcSaveSheetXVFromMapQ(Save_Sheetx, Save_Sheetv, m_mapQ);
        
        calcHitHistFftMap(hitFFTMap, hitmap);
    }
    
    //-----Q Threshold = 45ADC = 1.5fC = 1.5*1e-15 C --------
    //-----Q Satruation = 45fC
    //-----Q Resolution = 0.256fC
    //-----T jitter = 2ns
 
    double Q_threshold  = 1.5e-15 / 1.602176565e-19;
    double Q_saturation = 45.e-15 / 1.602176565e-19;
    double Q_resolution = 0.256e-15 / 1.602176565e-19;
    double T_threshold  = 12; // 12 mV
    //double T_threshold = 24; // 24 mV
    //double T_threshold = 18; // 18 mV
 
    // cout << "---------------------" << endl;
 
    //m000_nXstrips = m_mapQ[0][0].size() + m_mapQ[1][0].size();
    //m000_nVstrips = m_mapQ[0][1].size() + m_mapQ[1][1].size();
 
    //cout << m000_nXstrips << "  " << Nx_below_threshold << "  " << m000_nVstrips << "  " << Nv_below_threshold << endl;
    
    for (int i = 0; i < Save_Sheetx.size(); i++) {
        int sheetx_id = Save_Sheetx[i] / 10000;
        int stripx_id = Save_Sheetx[i] % 10000;
        //if(m_mapQ[sheetx_id][0][stripx_id]>=Q_threshold){
        Save_FinalstripIDx.push_back(Save_Sheetx[i]);
        //std::cout << "zhaojy check InductionGar2 sheet -> " << m000_xstripSheet.size() << " " << sheetx_id << "  " << stripx_id << std::endl;
        m000_xstripSheet.push_back(sheetx_id);
        m000_xstripID.push_back(stripx_id);
        m000_xstripQ.push_back(m_mapQ[sheetx_id][0][stripx_id]);
        //}
    }
 
    for (int i = 0; i < Save_Sheetv.size(); i++) {
        int sheetv_id = Save_Sheetv[i] / 10000;
        int stripv_id = Save_Sheetv[i] % 10000;
        //if(m_mapQ[sheetv_id][1][stripv_id]>=Q_threshold){
        Save_FinalstripIDv.push_back(Save_Sheetv[i]);
        m000_vstripSheet.push_back(sheetv_id);
        m000_vstripID.push_back(stripv_id);
        m000_vstripQ.push_back(m_mapQ[sheetv_id][1][stripv_id]);
        //}
    }
 
    double Tmin, Tmin2;
    std::vector<double> xTfirst; //time when first electron leave 3rd gem
    std::vector<double> vTfirst;
    //T_Branch & E_Branch
    for (int h = 0; h < Save_FinalstripIDx.size(); h++) {
 
        std::map<int, std::vector<int> > electron_hit_tmp;
        const int _low                    = 0;
        const int _up                     = 1000;
        const int _nbins                  = 2000;
        double binsize                    = (double)(_up - _low) / _nbins;
        int sheetx_id                     = Save_FinalstripIDx[h] / 10000;
        int stripx_id                     = Save_FinalstripIDx[h] % 10000;
        double histX_bin_Content[_nbins]  = {0};
        double histX_bin_Content2[_nbins] = {0};
        Tmin = Tmin2 = 999999;
        if (isDeadStrip(layer,sheetx_id,0,stripx_id)) {
            xTfirst.push_back(Tmin);
            double T_th = -99.;
            m000_xstripT_Branch.push_back(T_th);
            double Qin = -99.;
            m000_xstripQ_Branch.push_back(Qin);
            continue;
        };
        // we might want to change this into a easier version
        // input: Electrons, etc.........
        // out: the current of this strip
        if (bShowRef) {
            int original_affected_elec = 0;
 
            for (int i = 0; i < Save_Gridx.size(); i++) {
                //int z_grid_x = Save_Gridx[i]/100;
                //int z_grid_v = Save_Gridx[i]%100/10;
                //int z_index  = Save_Gridx[i]%10;
                int z_IDx     = Save_IDx[i];
                int start_bin = Save_Tx[i] / binsize;
                // std::cout << "start_bin  x" << start_bin << std::endl;
 
                if (z_IDx == Save_FinalstripIDx[h]) {
                    if (Save_Tx[i] < Tmin) Tmin = Save_Tx[i];
                    if (start_bin < _nbins and start_bin >= 0) {
                        electron_hit_tmp[Save_Gridx[i]].push_back(start_bin); //prepare the hit list. should have boundary check
                                                                              //discards evt >= _nbins or evt < 0.
                        original_affected_elec++;
                    }
                    //for(int addi = start_bin; addi < start_bin+160; addi+=2){
                    //  histX_bin_Content[addi]+=SignalX[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
                    //  histX_bin_Content[addi+1]+=SignalX[layer][z_grid_x][z_grid_v][z_index][(addi-start_bin)/2+1];
                    //}
                }
            }
 
            for (std::map<int, std::vector<int> >::iterator it = electron_hit_tmp.begin();
                 it != electron_hit_tmp.end(); it++) {
                std::vector<int> &hitlist = it->second;
                const int &Save_Grid      = it->first;
                // legacy
                if (hitlist.size() < electrons_select_method_threhold) {
                    conv1PerGrid_legacy(Save_Grid, histX_bin_Content, hitlist, layer, 0);
                } else {
                    conv1PerGrid_fft(Save_Grid, histX_bin_Content, hitlist, layer, 0);
                }
            }
 
            cout << "original_affected_elec " << original_affected_elec << "  ";
            calcStripCurrent(histX_bin_Content2, Tmin2,hitFFTMap, hitmap, layer, sheetx_id, stripx_id, 0);
 
            if (Tmin != Tmin2) {
                cout << "at stripX " << stripx_id << "Tmin: " << Tmin << " Tmin2: " << Tmin2 << "\n";
            }
            { // compare results 
 
                //histX_bin_Content2 and histX_bin_Content
                compareArray(histX_bin_Content, histX_bin_Content2, _nbins, "hist_bin_content", 1e-5);
                string emptyStr = "";
                DrawToFile(histX_bin_Content, (emptyStr + "st" + (long)stripx_id + "old").Data());
                DrawToFile(histX_bin_Content2, (emptyStr + "st" + (long)stripx_id + "new").Data());
            }
 
            // electronic parts
 
        } else {
            //calcStripCurrent(histX_bin_Content, Tmin, hitmap, layer, sheetx_id, stripx_id, 0);
            calcStripCurrent(histX_bin_Content2, Tmin,hitFFTMap, hitmap, layer, sheetx_id, stripx_id, 0);
            //compareArray(histX_bin_Content,histX_bin_Content2, _nbins,"hitHistFftMethod", 1e-5);
        }
 
        xTfirst.push_back(Tmin);
 
        #pragma region
        //TH1F h_signalT = Convolution_Tbranch(histX_bin_Content);
        TH1D h_signalT = ctf.Convolution_Tbranch_fft(histX_bin_Content);
 
        int flg_xT  = 0;
        double T_th = 0.;
        T_threshold = T_thr_V[layer][sheetx_id][0][stripx_id];
        //if(T_threshold<=0) T_threshold=thresholdMin_X_T;
        //T_threshold = 46.7;
        for (int init = 1; init < _nbins; init++) {
            if (flg_xT == 0 && fabs(h_signalT.GetBinContent(init + 1)) >= fabs(T_threshold)) {
                T_th = _low + binsize * (init - 1) + fabs(T_threshold - h_signalT.GetBinContent(init)) * binsize / fabs(h_signalT.GetBinContent(init + 1) - h_signalT.GetBinContent(init)) + gRandom->Gaus(0, 2);
                m000_xstripT_Branch.push_back(T_th); // jitter 2ns
                flg_xT = 1;
            }
        }
        if (flg_xT == 0) {
            T_th = -99.;
            m000_xstripT_Branch.push_back(T_th);
        }
        double V_sampling = 0.;
        double Qin        = 0.;
        double Efine      = 0.;
 
        if (T_th > 0) {
            //TH1D h_signalE = Convolution_Ebranch(histX_bin_Content);//debug; do not need to do if V_T(max) < T_threshold(mV)
            TH1D h_signalE     = cef.Convolution_Ebranch_fft(histX_bin_Content); //debug; do not need to do if V_T(max) < T_threshold(mV)
            double T_sample    = T_th + sample_delay;                            //ns
            int sample_bin     = T_sample / binsize;
            double sample_bin_ = T_sample / binsize;
            V_sampling         = h_signalE.GetBinContent(sample_bin) + (h_signalE.GetBinContent(sample_bin + 1) - h_signalE.GetBinContent(sample_bin)) * double(sample_bin_ - sample_bin);
            int over_th        = 0;
 
 
            if (h_signalE.GetMaximum() > E_thr_V[layer][sheetx_id][0][stripx_id]) over_th = 1;
 
            if (over_th == 1) {
                Qin = V_sampling*VQ_slope + VQ_const+gRandom->Gaus(0,m_QinGausSigma[0]);
                if (m_saturation && Qin > Qsaturation[layer][sheetx_id][0][stripx_id])
                    Qin = Qsaturation[layer][sheetx_id][0][stripx_id];
            } else
                Qin = -99;
        } else {
            V_sampling = -99.;
            Qin        = -99.;
            //m000_xstripQ_Branch.push_back(Qin);
        }
 
        m000_xstripQ_Branch.push_back(Qin);
        #pragma endregion
 
    }
 
    for (int h = 0; h < Save_FinalstripIDv.size(); h++) {
        std::map<int, std::vector<int> > electron_hit_tmp;
        const int _low                   = 0;
        const int _up                    = 1000;
        const int _nbins                 = 2000;
        double binsize                   = (double)(_up - _low) / _nbins;
        int sheetv_id                    = Save_FinalstripIDv[h] / 10000;
        int stripv_id                    = Save_FinalstripIDv[h] % 10000;
        double histV_bin_Content[_nbins] = {0};
        Tmin                             = 999999;
 
        if (isDeadStrip(layer,sheetv_id,1,stripv_id)) {
            vTfirst.push_back(Tmin);
            double T_th = -99.;
            m000_vstripT_Branch.push_back(T_th);
            double Qin = -99.;
            m000_vstripQ_Branch.push_back(Qin);
            continue;
        };
 
        calcStripCurrent(histV_bin_Content, Tmin, hitFFTMap,hitmap, layer, sheetv_id, stripv_id, 1);
 
        vTfirst.push_back(Tmin);
 
        TH1D h_signalT = ctf.Convolution_Tbranch_fft(histV_bin_Content);
 
        int flg_vT  = 0;
        double T_th = 0.;
        T_threshold = T_thr_V[layer][sheetv_id][1][stripv_id];
        //if(T_threshold<=0) T_threshold=thresholdMin_V_T;
        //T_threshold = 46.7;
        for (int init = 1; init < _nbins; init++) {
            if (flg_vT == 0 && fabs(h_signalT.GetBinContent(init + 1)) >= fabs(T_threshold)) {
                T_th = _low + binsize * (init - 1) + fabs(T_threshold - h_signalT.GetBinContent(init)) * binsize / fabs(h_signalT.GetBinContent(init + 1) - h_signalT.GetBinContent(init)) + gRandom->Gaus(0, 2);
                m000_vstripT_Branch.push_back(T_th); // jitter 2ns
                flg_vT = 1;
            }
        }
        if (flg_vT == 0) {
            T_th = -99.;
            m000_vstripT_Branch.push_back(T_th);
        }
        double V_sampling = 0.;
        double Qin        = 0.;
        double Efine      = 0.;
        if (T_th > 0) {
            //TH1F h_signalE = Convolution_Ebranch(histV_bin_Content);
            TH1D h_signalE     = cef.Convolution_Ebranch_fft(histV_bin_Content);
            double T_sample    = T_th + sample_delay; //ns
            int sample_bin     = T_sample / binsize;
            double sample_bin_ = T_sample / binsize;
            V_sampling         = h_signalE.GetBinContent(sample_bin) + (h_signalE.GetBinContent(sample_bin + 1) - h_signalE.GetBinContent(sample_bin)) * double(sample_bin_ - sample_bin);
            int over_th        = 0;
 
            if (h_signalE.GetMaximum() > E_thr_V[layer][sheetv_id][1][stripv_id]) over_th = 1;
 
            if (over_th == 1) {
                Qin = V_sampling*VQ_slope + VQ_const+gRandom->Gaus(0,m_QinGausSigma[1]);
                if (m_saturation && Qin > Qsaturation[layer][sheetv_id][1][stripv_id])
                    Qin = Qsaturation[layer][sheetv_id][1][stripv_id];
            } else
                Qin = -99;
        } else {
            V_sampling = -99.;
            Qin        = -99.;
            //m000_vstripQ_Branch.push_back(Qin);
        }
 
        m000_vstripQ_Branch.push_back(Qin);
        // cout<<"V    layer:"<<layer<<"    sheet:"<<sheetv_id<<"    strip:"<<stripv_id<<"    T_Branch:"<<T_th<<"    V_sample:"<<V_sampling<<endl;
        /*  
            int layer_tem = layer;
            char temp0[1]; sprintf(temp0, "%i", layer_tem);
            char temp1[1]; sprintf(temp1, "%i", sheetv_id);
            char temp2[3]; sprintf(temp2, "%i", stripv_id);
            string FILEname = FilePath + "/share/output/V-" + temp0 + '-' + temp1 + '-' + temp2 + ".root";
            TFile* f_x1  = new TFile(FILEname.c_str(), "recreate"); f_x1->cd(); sig_current.Write(); h_signalT.Write();  f_x1->Close(); delete f_x1;
            */
    }
 
    //-------------------------------------------------------------------------------
    //cout<<"push back"<<endl;
    double q_to_fC_factor = 1.602176565e-4;
    //std::cout << "loop V_Q finished " << std::endl;
 
    //cout<<"X"<<endl;
    if (storeFlag) {
        m_nXstrips = 0;
        for (int i = 0; i < m000_xstripSheet.size(); i++) {
            m_xstripSheet.push_back(m000_xstripSheet[i]);
            m_xstripID.push_back(m000_xstripID[i]);
            m_xstripQ.push_back(m000_xstripQ[i] * q_to_fC_factor);
            m_xstripT_Branch.push_back(m000_xstripT_Branch[i]);
            m_xstripQ_Branch.push_back(m000_xstripQ_Branch[i]);
            m_xTfirst.push_back(xTfirst[i]);
            m_nXstrips++;
        }
        m_nVstrips = 0;
        //cout<<"V"<<endl;
        for (int i = 0; i < m000_vstripSheet.size(); i++) {
            m_vstripSheet.push_back(m000_vstripSheet[i]);
            m_vstripID.push_back(m000_vstripID[i]);
            m_vstripQ.push_back(m000_vstripQ[i] * q_to_fC_factor);
            m_vstripT_Branch.push_back(m000_vstripT_Branch[i]);
            m_vstripQ_Branch.push_back(m000_vstripQ_Branch[i]);
            m_vTfirst.push_back(vTfirst[i]);
            m_nVstrips++;
        }
    } else {
        m_nXstrips = 0;
        for (int i = 0; i < m000_xstripSheet.size(); i++) {
            if (m000_xstripQ_Branch[i] >= 0) {
                m_xstripSheet.push_back(m000_xstripSheet[i]);
                m_xstripID.push_back(m000_xstripID[i]);
                m_xstripQ.push_back(m000_xstripQ[i] * q_to_fC_factor);
                m_xstripT_Branch.push_back(m000_xstripT_Branch[i]);
                m_xstripQ_Branch.push_back(m000_xstripQ_Branch[i]);
                m_xTfirst.push_back(xTfirst[i]);
                m_nXstrips++;
                //              cout<<m000_xstripQ[i]*q_to_fC_factor<<"    "<<m000_xstripQ_Branch[i]<<endl;
            }
        }
        m_nVstrips = 0;
        //cout<<"V"<<endl;
        for (int i = 0; i < m000_vstripSheet.size(); i++) {
            if (m000_vstripQ_Branch[i] >= 0) {
                m_vstripSheet.push_back(m000_vstripSheet[i]);
                m_vstripID.push_back(m000_vstripID[i]);
                m_vstripQ.push_back(m000_vstripQ[i] * q_to_fC_factor);
                m_vstripT_Branch.push_back(m000_vstripT_Branch[i]);
                m_vstripQ_Branch.push_back(m000_vstripQ_Branch[i]);
                m_vTfirst.push_back(vTfirst[i]);
                m_nVstrips++;
                //                  cout<<m000_vstripQ[i]*q_to_fC_factor<<"    "<<m000_vstripQ_Branch[i]<<endl;
            }
        }
    }
    //checkMemorySize();
 
    //cout<<"InductionGar2::setMultiElectrons() ends"<<endl;
}
 
/**
 * @brief calcs m_mapQ, Q integrated from
 * 
 * @param histmap 
 * @param layer 
 */
void InductionGar2::calcMapQFromHitHistMap(const HitHistMap &histmap, int layer, std::map<int, double> mapQ[2][2]) const {
    for (HitHistMap::const_iterator it = histmap.begin(); it != histmap.end(); ++it) {
        const IndexGar &index = it->first;
        const Vint &hist      = it->second;
        int xID               = index.stripX;
        int vID               = index.stripV;
        int grid_x            = index.gridX;
        int grid_v            = index.gridV;
        int sheet             = index.sheet;
        if (xID < 0 || vID < 0) continue;
        CgemGeoReadoutPlane *CgemReadoutPlane = m_geomSvc->getReadoutPlane(layer, sheet);
        int nXID                              = CgemReadoutPlane->getNXstrips();
        int nVID                              = CgemReadoutPlane->getNVstrips();
        int nElectrons                        = 0;
        for (int i = 0; i < _nbins + 1; i++) { // last bin mean overflow.
            nElectrons += hist[i];
        }
        for (int indexOID = 0; indexOID < 4; indexOID++) { // X
            int dID = indexOID - 2;                        //i.e. indexOID==2 show it is on this strip
            if (0 <= xID + dID && xID + dID < nXID) {
                mapQ[sheet][0][xID + dID] += RatioX[layer][grid_x][grid_v][indexOID] * nElectrons;
            }
        }
        for (int indexOID = 0; indexOID < 3; indexOID++) { // V
            int dID = indexOID - 1;                        //i.e. indexOID==1 show it is on this strip
            if (0 <= vID + dID && vID + dID < nVID) {
                mapQ[sheet][1][vID + dID] += RatioV[layer][grid_x][grid_v][indexOID] * nElectrons;
            }
        }
    } // may contain Q==0.
    // as far as i know, numbering scheme for XV strip is like {{1,...,n_max_sheet1},{n_max_sheet1+1,...n_max_sheet2}} (CITE:https://docbes3.ihep.ac.cn/~offlinesoftware/images/d/da/CGEM-IT_geometry_from_the_final_design.pdf)
    //is but here, and in Induction and GeomSvc, the ID of strips appears to begin from 0 in its both sheet. OK.
}
void InductionGar2::calcSaveSheetXVFromMapQ(Vint &Save_Sheetx, Vint &Save_Sheetv, std::map<int, double> const mapQ[2][2]) const {
    //std::map<int, double> m_mapQ[2][2]; //[nSheet][XV-view] {index of ID, Q}
    for (int sheet = 0; sheet < 2; sheet++) {
        for (int xvView = 0; xvView < 2; xvView++) {
            Vint &save_sheet = xvView ? Save_Sheetv : Save_Sheetx; //0 as x and 1 as v
            for (std::map<int, double>::const_iterator it = mapQ[sheet][xvView].begin();
                 it != mapQ[sheet][xvView].end(); it++) {
                if (it->second != 0) {
                    save_sheet.push_back(sheet * 10000 + it->first);
                }
            }
        }
    }
}
/**
 * @brief calcs
 * 
 * @param strip_current the current of strip. this func will not initialize the current.
 * @param histmap 
 * @param layer 
 * @param req_sheetid 
 * @param req_stripid 
 * @param xvaxis 
 */
void InductionGar2::calcStripCurrent(double *strip_current, double &timeFirstHit, const HitHistMap &histmap, int layer, int req_sheetid, int req_stripid, int xvaxis) const {
 
    const int _low                                   = 0;
    const int _up                                    = 1000;
    const int _nbins                                 = 2000;
    double binsize                                   = (double)(_up - _low) / _nbins;
    ConvolveWithConst(&Signal_Convolver)[3][5][5][4] = const_cast<ConvolveWithConst(&)[3][5][5][4]>(xvaxis ? Signal_ConvolverV : Signal_ConvolverX);
    // x[ ][ ][m][ ]
    // v[ ][m][ ][0]
    int indexTimeFirstHit = 9999;
    int affectedElectron  = 0;
    for (HitHistMap::const_iterator it = histmap.begin(); it != histmap.end(); ++it) {
        const IndexGar &index = it->first;
        const Vint &hist      = it->second;
 
        int xID    = index.stripX;
        int vID    = index.stripV;
        int strip  = xvaxis ? vID : xID;
        int grid_x = index.gridX;
        int grid_v = index.gridV;
        int sheet  = index.sheet;
        int indexOID;
        if (xID < 0 || vID < 0) continue;
        //if (req_sheetid!=sheet || req_stripid != strip) continue;
        if (req_sheetid != sheet) continue;
        if (xvaxis == 0) { //x
            //indexOID=strip-req_stripid+2;//
            indexOID = req_stripid - strip + 2; //
            if (indexOID < 0 || indexOID >= 4) continue;
            if (RatioX[layer][grid_x][grid_v][indexOID] == 0) continue;
 
        } else { //v
            //indexOID=strip-req_stripid+1;
            indexOID = req_stripid - strip + 2;
            if (indexOID < 0 || indexOID >= 3) continue;
            if (RatioV[layer][grid_x][grid_v][indexOID] == 0) continue;
        }
        // now we get valid coordinates that interfere with the strip.
 
        double electron_hit_tmp_hist[_nbins] = {0};
        for (int i = 0; i < _nbins; i++) electron_hit_tmp_hist[i] = hist[i];
 
        double strip_current_tmp[_nbins] = {0};
        Signal_Convolver[layer][grid_x][grid_v][indexOID].convolve(strip_current_tmp, electron_hit_tmp_hist, 0, _nbins);
        for (int i = 0; i < _nbins; i++) {
            strip_current[i] += strip_current_tmp[i];
        }
 
        // for Tmin stuff.
        int indexTimeFirstHitLocal = 9999;
        for (int i = 0; i < _nbins; i++) {
            if (hist[i] != 0) {
                indexTimeFirstHitLocal = i;
                break;
            }
        }
        if (indexTimeFirstHitLocal < indexTimeFirstHit) indexTimeFirstHit = indexTimeFirstHitLocal;
 
        // affected Electrons
        for (int i = 0; i < _nbins; i++) {
            affectedElectron += hist[i];
        }
    }
    timeFirstHit = binsize * (indexTimeFirstHit); // const error: binsize
    //cout<<"laterFffectedElectron "<< affectedElectron<<endl;
}
void InductionGar2::calcStripCurrent(double *strip_current, double &timeFirstHit, const std::map<IndexGar,std::vector<double> > &histFFTmap,const HitHistMap &histmap, int layer, int req_sheetid, int req_stripid, int xvaxis) const {
 
    const int _low                                   = 0;
    const int _up                                    = 1000;
    const int _nbins                                 = 2000;
    double binsize                                   = (double)(_up - _low) / _nbins;
    const ConvolveWithConst(&Signal_Convolver)[3][5][5][4] = xvaxis ? Signal_ConvolverV : Signal_ConvolverX;
    // x[ ][ ][m][ ]
    // v[ ][m][ ][0]
    int indexTimeFirstHit = 9999;
    int affectedElectron  = 0;
 
    const int length_Inner=Signal_ConvolverX[0][0][0][0].getLength();
    const int HalfPlus1=(length_Inner/2+1);
    vector<double> stripCurrentFFT ;
    stripCurrentFFT.resize(HalfPlus1*2);
    double* pStripCurrentFFT=&stripCurrentFFT.front();
    for (HitHistMap::const_iterator it = histmap.begin(); it != histmap.end(); ++it) {
        const IndexGar &index = it->first;
        const Vint &hist      = it->second;
        
        const std::vector<double>  &histFft      =histFFTmap.at(index);
        const double * pHistFft=&histFft.front();
 
        int xID    = index.stripX;
        int vID    = index.stripV;
        int strip  = xvaxis ? vID : xID;
        int grid_x = index.gridX;
        int grid_v = index.gridV;
        int sheet  = index.sheet;
        int indexOID;
        if (xID < 0 || vID < 0) continue;
        //if (req_sheetid!=sheet || req_stripid != strip) continue;
        if (req_sheetid != sheet) continue;
        if (xvaxis == 0) { //x
            //indexOID=strip-req_stripid+2;//
            indexOID = req_stripid - strip + 2; //
            if (indexOID < 0 || indexOID >= 4) continue;
            if (RatioX[layer][grid_x][grid_v][indexOID] == 0) continue;
 
        } else { //v
            //indexOID=strip-req_stripid+1;
            indexOID = req_stripid - strip + 2;
            if (indexOID < 0 || indexOID >= 3) continue;
            if (RatioV[layer][grid_x][grid_v][indexOID] == 0) continue;
        }
        // now we get valid coordinates that interfere with the strip.
 
 
        ////double electron_hit_tmp_hist[_nbins] = {0};
        ////for (int i = 0; i < _nbins; i++) electron_hit_tmp_hist[i] = hist[i];
//      
        Signal_Convolver[layer][grid_x][grid_v][indexOID].MultiplyAndAdd((ConvolveWithConst::Complex*)pStripCurrentFFT,(ConvolveWithConst::Complex*)pHistFft);
        //if (compareArray(pStripCurrentFFT,zeros,_nbins/2,"",1e-10,true))cout<<"CgemDigi...::Induc...r2::calcCurrent: ffound all bin empty!!!! "<<endl;
        //double strip_current_tmp[_nbins] = {0};
        //Signal_Convolver[layer][grid_x][grid_v][indexOID].convolve(strip_current_tmp, electron_hit_tmp_hist, 0, _nbins);
        //for (int i = 0; i < _nbins; i++) {
        //    strip_current[i] += strip_current_tmp[i];
        //}
 
 
        // for Tmin stuff.
        int indexTimeFirstHitLocal = 9999;
        for (int i = 0; i < _nbins; i++) {
            if (hist[i] != 0) {
                indexTimeFirstHitLocal = i;
                break;
            }
        }
        if (indexTimeFirstHitLocal < indexTimeFirstHit) indexTimeFirstHit = indexTimeFirstHitLocal;
 
        // affected Electrons
        for (int i = 0; i < _nbins; i++) {
            affectedElectron += hist[i];
        }
    }
    //if (compareArray(pStripCurrentFFT,zeros,_nbins/2,"",1e-10,true))cout<<"CgemDigi...::Induc...r2::calcCurrent aftersum: found all bin empty!! "<<endl;
    Signal_ConvolverX[0][0][0][0].IFFT(strip_current,(ConvolveWithConst::Complex*)pStripCurrentFFT,_nbins,HalfPlus1*2);
    //if (compareArray(strip_current,zeros,_nbins/2,"",1e-10,true))cout<<"CgemDigi...::Induc...r2::calcCurrent strip_current: found all bin empty!! "<<endl;
    timeFirstHit = binsize * (indexTimeFirstHit); // const error: binsize
    //cout<<"laterFffectedElectron "<< affectedElectron<<endl;
}
 
void InductionGar2::calcHitHistFftMap(std::map<IndexGar,std::vector<double> > &histFFTMap,const HitHistMap &histmap)const{
    
    //TVirtualFFT* fft=Signal_ConvolverX[0][0][0][0].getFFT();
    const int length_Inner=Signal_ConvolverX[0][0][0][0].getLength();
    vector<double> vIn;
    vIn.resize(_nbins);
    double* pIn=&vIn.front();
    //vector<double> vOut;
    //vOut.resize(length_Inner+2);
    //ConvolveWithConst::Complex* pcOut=(ConvolveWithConst::Complex*) &vOut.front();
    for (HitHistMap::const_iterator it = histmap.begin(); it != histmap.end(); ++it){
        IndexGar index=it->first;
        const Vint &hist=it->second;
        histFFTMap[index] .resize((length_Inner/2+1)*2);
        double * pOut=&histFFTMap[index].front();
        //int * pIn=&histmap[index].front();
        for(int i=0; i<_nbins; i++){
            vIn[i]=hist[i];
        }
        Signal_ConvolverX[0][0][0][0].FFT((ConvolveWithConst::Complex*)pOut,pIn,(length_Inner/2+1)*2,_nbins);
        //if (compareArray(pOut,zeros,_nbins/2,"",1e-10,true))cout<<"CgemDigi...::Induc...r2::calcHitHistFFT : found all bins empty at strip= "<<index.stripX<<": "<<index.stripV <<endl;
        //vIn.clear();
        //for (int i=0; i<_nbins; ++i){
        //    vIn[i]=it->second[i];
        //}
        //fft->SetPoints(pIn);
        //fft->Transform();
        //fft->GetPoints(pOut);
    }
}
//i may want to warp this into a function
/**
 * @brief legacy method of doing conv1 at a single grid for a strip. additive on hist_bin_Content
 *  discards t > _nbins or t <0  events
 * @param Save_Grid [grid_x]*100+[grid_v]*10+[z_index]*1
 * @param hist_bin_content to be summed on
 * @param hitlist list of hits;
 * @param layer layer.
 * @param axis 0:X, 1:V
 */
void InductionGar2::conv1PerGrid_legacy(int Save_Grid, double *hist_bin_Content,
                                        const std::vector<int> &hitlist, int layer, int axis) const {
    //double ***** Signal=0;
    // this should not work, since you get no dimention info
    // and this following one might work.
 
    double const(&Signal)[3][5][5][4][100] = axis ? SignalV : SignalX;
    // our old C++ compiler cling(from old ROOT) complains about this.
    // thankfully, our g++ does not.
 
    const int _nbins = 2000;
    int z_grid_x     = Save_Grid / 100;
    int z_grid_v     = Save_Grid % 100 / 10;
    int z_index      = Save_Grid % 10;
    for (std::vector<int>::const_iterator it_start_bin = hitlist.begin(); it_start_bin != hitlist.end(); it_start_bin++) {
        if (*it_start_bin < 0) continue;
        for (int addi = *it_start_bin; (addi < *it_start_bin + 160) and (addi < _nbins - 1); addi += 2) {         // 160 cell-layer move. i believe more hits in a semi-cell make fft-convolve relatively faster.
            hist_bin_Content[addi] += Signal[layer][z_grid_x][z_grid_v][z_index][(addi - *it_start_bin) / 2 + 1]; // caution, its indice starts from...1?
            hist_bin_Content[addi + 1] += Signal[layer][z_grid_x][z_grid_v][z_index][(addi - *it_start_bin) / 2 + 1];
        } // here its content is like stage with step size 2;
    }
}
 
/**
 * @brief fft method of doing conv1 at a single grid for a strip. same usage as conv1PerGrid_legacy
 * discards t > _nbins or t <0  events
 * @param Save_Grid [grid_x]*100+[grid_v]*10+[z_index]*1
 * @param hist_bin_content to be summed on
 * @param hitlist list of hits
 * @param layer layer.
 * @param axis 0:X, 1:V
 */
void InductionGar2::conv1PerGrid_fft(int Save_Grid, double *hist_bin_Content,
                                     const std::vector<int> &hitlist, int layer, int axis) const {
    ConvolveWithConst(&Signal_Convolver)[3][5][5][4] = const_cast<ConvolveWithConst(&)[3][5][5][4]>(axis ? Signal_ConvolverV : Signal_ConvolverX);
 
    const int _nbins        = 2000;
    double hist_tmp[_nbins] = {0};
    int z_grid_x            = Save_Grid / 100;
    int z_grid_v            = Save_Grid % 100 / 10;
    int z_index             = Save_Grid % 10;
 
    //bool found_mismatch          = false;
    double electron_hit_tmp_hist[_nbins] = {0};
 
    for (std::vector<int>::const_iterator it_start_bin = hitlist.begin(); it_start_bin != hitlist.end(); it_start_bin++) {
        if (*it_start_bin < _nbins and *it_start_bin >= 0) electron_hit_tmp_hist[*it_start_bin] += 1;
    }
    Signal_Convolver[layer][z_grid_x][z_grid_v][z_index].convolve(hist_tmp, electron_hit_tmp_hist, 0, _nbins);
 
#ifdef INDUCTIONGAR_DEBUG_CONV1PERGRID_FFT_TEST
    {
        double hist_tmp2[_nbins]          = {0};
        bool found_mismatch               = false;
        const double Accept_Threshold     = 1e-5;
        const double Accept_Threshold_ref = 1e-5;
        conv1PerGrid_legacy(Save_Grid, hist_tmp2, hitlist, layer, axis);
        for (int i = 0; i < _nbins; i++) {
            if (abs(hist_tmp[i] - hist_tmp2[i]) >= Accept_Threshold and abs((hist_tmp[i] - hist_tmp2[i]) / hist_tmp2[i]) >= Accept_Threshold_ref) {
                if (not found_mismatch) {
                    std::cout << "CgemDigitizerSvc::InductionGar2::conv1PerGrid_fft found results not match:  " << endl;
                    found_mismatch = true;
                }
                std::cout << "    previous: " << hist_tmp[i] << "; new " << hist_tmp2[i] << "; at i=" << i << '\n'
            }
        }
    }
#endif
 
    //add
    for (int i = 0; i < _nbins; i++) {
        hist_bin_Content[i] += hist_tmp[i];
    }
}
 
/**
 * @brief fft method of doing conv1 at a single grid for a strip.
 * 
 * @param gridX 
 * @param gridV 
 * @param indexOID 
 * @param hist_bin_Content 
 * @param elecHitHist 
 * @param layer 
 * @param axis 
 */
void InductionGar2::conv1PerGrid_legacy(int gridX, int gridV, int indexOID, double *hist_bin_Content,
                                        const std::vector<int> &elecHitHist, int layer, int axis) const {
    //double ***** Signal=0;
    // this should not work, since you get no dimention info
    // and this following one might work.
 
    const double(&Signal)[3][5][5][4][100] = axis ? SignalV : SignalX;
 
    // our old C++ compiler cling(from old ROOT) complains about this.
    // thankfully, our g++ does not.
 
    const int _nbins = 2000;
 
    for (int start_bin = 0; start_bin < _nbins; start_bin++) {
        int nElec = elecHitHist[start_bin];
        for (int addi = start_bin; (addi < start_bin + 160) and (addi < _nbins - 1); addi += 2) {
            hist_bin_Content[addi] += Signal[layer][gridX][gridV][indexOID][(addi - start_bin) / 2 + 1] * nElec;
            hist_bin_Content[addi + 1] += Signal[layer][gridX][gridV][indexOID][(addi - start_bin) / 2 + 1] * nElec;
        } // a step2 thing.
    }
 
    //int z_grid_x     = Save_Grid / 100;
    //int z_grid_v     = Save_Grid % 100 / 10;
    //int z_index      = Save_Grid % 10;
    //for (std::vector<int>::const_iterator it_start_bin = hitlist.begin(); it_start_bin != hitlist.end(); it_start_bin++) {
    //    if (*it_start_bin < 0) continue;
    //    for (int addi = *it_start_bin; (addi < *it_start_bin + 160) and (addi < _nbins - 1); addi += 2) {         // 160 cell-layer move. i believe more hits in a semi-cell make fft-convolve relatively faster.
    //        hist_bin_Content[addi] += Signal[layer][z_grid_x][z_grid_v][z_index][(addi - *it_start_bin) / 2 + 1]; // caution, its indice starts from...1?
    //        hist_bin_Content[addi + 1] += Signal[layer][z_grid_x][z_grid_v][z_index][(addi - *it_start_bin) / 2 + 1];
    //    } // here its content is like stage with step size 2;
    //}
}
 
void InductionGar2::conv1PerGrid_fft(int gridX, int gridV, int indexOID, double *hist_bin_Content,
                                     const std::vector<int> &elecHitHist, int layer, int axis) const {
    ConvolveWithConst(&Signal_Convolver)[3][5][5][4] = const_cast<ConvolveWithConst(&)[3][5][5][4]>(axis ? Signal_ConvolverV : Signal_ConvolverX);
 
    const int _nbins        = 2000;
    double hist_tmp[_nbins] = {0};
    int z_grid_x            = gridX;
    int z_grid_v            = gridV;
    int z_index             = indexOID;
 
    //bool found_mismatch          = false;
    double electron_hit_tmp_hist[_nbins] = {0};
    for (int i = 0; i < _nbins; i++) {
        electron_hit_tmp_hist[i] = elecHitHist[i];
    }
    //
    //for (std::vector<int>::const_iterator it_start_bin = hitlist.begin(); it_start_bin != hitlist.end(); it_start_bin++) {
    //    if (*it_start_bin < _nbins and *it_start_bin >= 0) electron_hit_tmp_hist[*it_start_bin] += 1;
    //}
 
    Signal_Convolver[layer][z_grid_x][z_grid_v][z_index].convolve(hist_tmp, electron_hit_tmp_hist, 0, _nbins);
 
#ifdef INDUCTIONGAR_DEBUG_CONV1PERGRID_FFT_TEST
    {
        double hist_tmp2[_nbins]          = {0};
        bool found_mismatch               = false;
        const double Accept_Threshold     = 1e-10;
        const double Accept_Threshold_ref = 1e-5;
        conv1PerGrid_legacy(gridX, gridV, indexOID, hist_tmp2, elecHitHist, layer, axis);
        for (int i = 0; i < _nbins; i++) {
            if (abs(hist_tmp[i] - hist_tmp2[i]) >= Accept_Threshold and abs((hist_tmp[i] - hist_tmp2[i]) / hist_tmp2[i]) >= Accept_Threshold_ref) {
                if (not found_mismatch) {
                    std::cout << "CgemDigitizerSvc::InductionGar2::conv1PerGrid_fft found results not match:  " << endl;
                    found_mismatch = true;
                }
                std::cout << "    previous: " << hist_tmp[i] << "; new " << hist_tmp2[i] << "; at i=" << i << '\n'
            }
        }
    }
#endif
 
    //add
    for (int i = 0; i < _nbins; i++) {
        hist_bin_Content[i] += hist_tmp[i];
    }
}
 
 
 
/**
 * @brief checks if all entries in mapQ1 are same in mapQ2.
 * 
 * @param mapQ1 
 * @param mapQ2 
 */
static void compareMap(const std::map<int, double> mapQ1[2][2], const std::map<int, double> mapQ2[2][2]) {
 
    for (int sheet = 0; sheet < 2; sheet++) {
        for (int xvView = 0; xvView < 2; xvView++) {
            const std::map<int, double> &m1 = mapQ1[sheet][xvView];
            const std::map<int, double> &m2 = mapQ2[sheet][xvView];
            cout << "CgemDigitizerSvc::InductionGar2::compareMap current sheet " << sheet << "; xvView " << xvView << endl;
            if (m1.size() != m2.size()) {
                cout << "\tcompareMap: size of maps not match: map1: " << m1.size() << ";\tmap2: " << m2.size() << '\n';
            }
            for (std::map<int, double>::const_iterator it = m2.begin(); it != m2.end(); it++) {
                int key      = it->first;
                double value = it->second;
                try {
                    if (value != m1.at(key)) {
                        std::cout << "\tcompareMap: value at key " << key << " not match: map1:" << m1.at(key) << "\tmap2:" << m2.at(key) << '\n';
                    }
                } catch (const std::out_of_range &e) {
                    std::cout << "\tcompareMap: key " << key << "in map2 " << m2.at(key) << " has no coresponding entry in map1 \n";
                }
            }
        }
    }
}
 
static void DrawToFile(const double *h, const char *filename) {
    TH1D h1("h_InductionGar2", "", _nbins, 0, 1000);
    for (int i = 0; i < _nbins; i++) {
        h1.SetBinContent(i + 1, h[i]);
    }
    TCanvas cx1("cx1", "", 1);
    h1.Draw();
 
    string filebase = "save/";
    string fn       = filebase + filename + ".png";
    cx1.Print(fn.c_str());
}
 
 
 
static void checkMemorySize()
{
    ProcInfo_t info;
    gSystem->GetProcInfo(&info);
    double currentMemory = info.fMemResident/1024;// MB
    cout<<"CgemDigitizerSvc::InductionGar2: current memory size "<<currentMemory<<" MB"<<endl;
}
 
 
bool InductionGar2::isDeadStrip(int layer, int sheet, int XVview, int strip)
{
    std::set<int> &deadStrips= m_deadStrip[layer][sheet][XVview];
    //cout<<"layer  "<<layer<<"    sheet   "<<sheet<<"  XVview  "<<XVview<<"  strip  "<<strip;
    if (deadStrips.count(strip)==1){
        //cout<<"is dead strip"<<endl;
        return true;
    }
    //cout<<"is not dead"<<endl;
    return false;
}