Dynamic NFOUR option and clean up several classes

[MyoungchulK]

Dear reviewers, This PR is not just including update, but also has a lots of code cleaning. It can be overwhelming (possibly frustrating) to check all of them. I try to explain most of them in this document but I realized document itself also become long. I will be really appreciate, if you spare little bit of patient for this PR (and for me). Thank you!:)

1. Dynamic NFOUR option.

After we generate an E-field WF and apply a detector response, we interpolate the signal with pre/user-defined time width (TIMESTEP) and zero-pad the signal by also pre/user-defined pad length (NFOUR). Then, we moved on to the noise generation part and convolved them together.

If the signal strength is enormous or the distance to the antenna cluster is too close, unfortunately, the pad length is too small to contain the entire signal. And signal WF from the root output has an unnatural or clipped shape. We traditionally set the pad length (NFOUR) to be large (1024 to 2048) based on experience, but I found, It can still cause the clipped WF.

In this update, I added the DYNAMIC_NFOUR option in the Setting class. If the DYNAMIC_NFOUR is set to 1, the pad length will be automatically adjusted based on the original signal WF length. And prevent clipped WF. I must tell you this DYNAMIC_NFOUR option is only implemented in the time-domain mode (SIMULATION_MODE = 1). If user select the DYNAMIC_NFOUR mode without selecting time-domain mode, It will automatically use default NFOUR value.

DYNAMIC_NFOUR option is placed just after sim apply detector response.

/*! 
    MK added -2023-05-20-
    If user selected DYNAMIC_NFOUR to 1,
    Changes the length of zero-pad (NFOUR / 2) to have a same length with original time-domain signal
    These automatic changes will prevent the cliped WF issue no matter how signal is long
    User don't have to empirically adjust the NFOUR values
*/
int new_NFOUR = settings1->NFOUR; ///< Default is user/pre-configured NFOUR
int new_init_T = init_T;
if (settings1->DYNAMIC_NFOUR == 1) {
    new_NFOUR = pad_len; ///< same with original signal WF
    if (new_NFOUR < settings1->NFOUR) new_NFOUR = settings1->NFOUR; ///< if length of signal is smaller than default NFOUR, just use default NFOUR
    new_init_T = settings1->TIMESTEP * -1.e9 * ((double) new_NFOUR / 4 + RandomTshift); ///< recalculate initial time
    //! We save the longest/smallest NFOUR for placing signal into noise pad at later.
    //! Use maximum NFOUR will prevent entire signal to be placed outside of noise pad
    //! Use minimum NFOUR will allow trigger scan to start the scan from the begining of all the signals in the noise pad
    if (new_NFOUR < new_NFOUR_min) new_NFOUR_min = new_NFOUR; 
    if (new_NFOUR > new_NFOUR_max) new_NFOUR_max = new_NFOUR;
}
//! Filled interpolated signal WF into zero-pad
double new_volts_forint[new_NFOUR / 2];
double new_T_forint[new_NFOUR / 2];
for(int n = 0; n < new_NFOUR / 2; n++){
    new_T_forint[n] = new_init_T + (double)n * settings1->TIMESTEP * 1.e9;
}
stations[i].strings[j].antennas[k].New_NFOUR[ray_sol_cnt] = new_NFOUR; ///< save in Report branch

//! do linear interpolation
//! changed to sinc interpolation Dec 2020 by BAC
Tools::SincInterpolation(pad_len, T_forfft, V_forfft, new_NFOUR / 2, new_T_forint, new_volts_forint);

new_NFOUR will set new zero-pad for signal WF (new_volts_forint, new_T_forint), and the interpolationed WF from sine interpolation will be stored in the zero-pad array
Different NFOUR values set by each signal WF will be saved in Report branch.

While setting NFOUR, we save the longest/smallest NFOUR for placing signal into noise pad at later.
We use maximum NFOUR (new_NFOUR_max) to prevent entire signal from being placed outside of the noise pad
And We use minimum NFOUR (new_NFOUR_min) to allow the trigger scan to start the scan from the beginning of all the signals in the noise pad

The maximum NFOUR (new_NFOUR_max) will be used to set max_total_bin to contain the longest signal in noise pad

//! calculate total number of bins we need to do trigger check
//! set the max_total_bin length by using new_NFOUR_max value to contain the longest signal in noise pad. MK added -2023-05-23-
max_total_bin = (stations[i].max_arrival_time - stations[i].min_arrival_time) / settings1->TIMESTEP + new_NFOUR_max * 3 + trigger->maxt_diode_bin;    // make more time

The minimum NFOUR (new_NFOUR_min) will be used to allow trigger scan to start the scan from the beginning of all the signals in the noise pad

//! Use new_NFOUR_min will allow trigger scan to start the scan from the begining of all the signals in the noise pad. MK added -2023-05-23-
trig_search_init = trigger->maxt_diode_bin + new_NFOUR_min;  // give some time shift for mimicing force trig events
trig_i = trig_search_init;

// save trig search bin info default (if no global trig, this value saved)
stations[i].total_trig_search_bin = max_total_bin - trig_search_init;

Because of different signal pad length for each solution, whether two ray solutions are connected or not will be identify by sum of half length of each pad

//! signal_dbin should be smaller than sum of half length of each signal WF to have connceted signals
int new_NFOUR_d = stations[i].strings[j].antennas[k].New_NFOUR[m - 1] / 4 + stations[i].strings[j].antennas[k].New_NFOUR[m] / 4;
if (signal_dbin[m - 1] < new_NFOUR_d)
{
    // if two ray_sol time delay is smaller than time window
    connect_signals.push_back(1);
    // cout<<"need two signal connection!!"<<endl;
}

The signal + noise convolution part (Select_Wave_Convlv_Exchange() 1st func. 2nd func, 3rd func) is modified to accept the different NFOUR values.
Since each signal WF length is different, we need to set the length of convolution pad by checking all the input signals length. We also need to consider the case that one of the signal range can cover other signals range. So, setting convolution pad like below is necessary.

void Report::Select_Wave_Convlv_Exchange(Settings *settings1, Trigger *trigger, Detector *detector, int signalbin1, int signalbin2, vector <double> &V1, vector <double> &V2, int *noise_ID, int ID, int StationIndex, int new_NFOUR1, int new_NFOUR2) {

    int bin_value;
    //! MK added -2023-05-23-
    //! set the edge of convolution pad by checking all the input signals length
    //! We do this in case second signal is covering entire first signal or vise versa
    int min_bin1 = signalbin1 - new_NFOUR1 / 4; ///< 1st signal
    int min_bin2 = signalbin2 - new_NFOUR2 / 4; ///< 2nd signal
    int min_bin = min_bin1;
    if (min_bin1 > min_bin2) min_bin = min_bin2; ///< if 2nd signal edge is smaller than 1st signal, we use 2nd signal edge for pad edge
    int max_bin1 = signalbin1 + new_NFOUR1 / 4;
    int max_bin2 = signalbin2 + new_NFOUR2 / 4;
    int max_bin = max_bin2;
    if (max_bin1 > max_bin2) max_bin = max_bin1; ///< if 1st signal edge is bigger than 2nd signal, we use 1st signal edge for pad edge
    int BINSIZE = max_bin - min_bin + 1;
    if (BINSIZE > int(settings1->NFOUR / 2) && BINSIZE <= settings1->NFOUR) BINSIZE = settings1->NFOUR; ///< Let just use double length. It is easier to apply diode
    double V_tmp[BINSIZE];
    for(int bin_tmp=0; bin_tmp<BINSIZE; bin_tmp++) {
        V_tmp[bin_tmp] = 0.;
    }

At last, because of the possibility of dynamic convolution pad length, I decided to create new function in Detector class that re-padding diode values based on convolution pad length. If convolution pad is not the same as default NFOUR value, It will re-pad the diode value and transform to frequency spectrum. And the convolution functions will use that array.

//! do myconvlv and replace the diode response array
//! Since signal WF length can be different by DYNAMIC_NFOUR option, based on the signal WF length we use different length of f-domain diode array. MK added -2023-05-23-
if (BINSIZE == int(settings1->NFOUR / 2)) trigger->myconvlv( V_total_forconvlv, BINSIZE, detector->fdiode_real, V_total_forconvlv); ///< use default diode
else if (BINSIZE == settings1->NFOUR) trigger->myconvlv( V_total_forconvlv, BINSIZE, detector->fdiode_real_double, V_total_forconvlv); ///< use default double length diode
else { ///< signal length is different than any default diode, we repad the diode and use it
        detector->get_NewDynamicDiodeModel(BINSIZE); ///< repad diode value
        trigger->myconvlv( V_total_forconvlv, BINSIZE, detector->fdiode_real_dynamic_databin, V_total_forconvlv);
    }

2. clean up

I clean up the Detector, Trigger, and Report classes. Unfortunately, it was a bit hard to implement the Dynamic NFOUR without cleaning.

In the AraSim class, I delete two array. They are taking space and do nothing.

double x_V[settings1->NFOUR/2];
double y_V[settings1->NFOUR/2];

In the Detector class, I mainly removed arrays that have names, including NFOUR.
It is not used anymore, but somehow sim still generates the value, saves it in the below arrays, and does nothing.
Removed arrays in the header file

vector <double> FilterGain_NFOUR;
vector <double> PreampGain_NFOUR;
vector <double> FOAMGain_NFOUR;
double GetFilterGain_NFOUR(int bin) { return FilterGain_NFOUR[bin]; }
double GetPreampGain_NFOUR(int bin) { return PreampGain_NFOUR[bin]; }
double GetFOAMGain_NFOUR(int bin) { return FOAMGain_NFOUR[bin]; }

I also removed some script lines that related to the above lines. Below lines are example

    // for NFOUR/2 t domain array
    double xfreq_NFOUR[settings1->NFOUR/4+1];   // array for FFT freq bin
    double ygain_NFOUR[settings1->NFOUR/4+1];   // array for gain in FFT bin

    df_fft = 1./ ( (double)(settings1->NFOUR/2) * settings1->TIMESTEP );

    for (int i=0;i<settings1->NFOUR/4+1;i++) {    // this one is for DATA_BIN_SIZE
        xfreq_NFOUR[i] = (double)i * df_fft / (1.E6); // from Hz to MHz
    }

    Tools::SimpleLinearInterpolation( N, xfreq, ygain, settings1->NFOUR/4+1, xfreq_NFOUR, ygain_NFOUR );

    for (int i=0;i<settings1->NFOUR/4+1;i++) {
        FilterGain_NFOUR.push_back( ygain_NFOUR[i] );
    }

I also added one more array and function in the detector class. When we start the simulation, we generate a diode response based on the NOUR length. In this update, I propose a dynamic NOUR option, in which NFOUR length will be changed by signal length (especially 1st and 2nd connected signal). So, when we convolved the noise and signal, we needed to re-pad the diode too. Below function will take the diode response from the sim and will re-pad when the signal length is different than NFOUR length.
Array in the header file
vector <double> fdiode_real_dynamic_databin;
Function

//! repadding diode value based on input pad length
void Detector::get_NewDynamicDiodeModel(int pad_len) {

    double diode_real_fft[pad_len * 2]; ///< double sized array for myconvlv

    for (int i = 0; i < pad_len * 2; i++) {  // 512 bin added for zero padding
        if ( i < (int)(maxt_diode / TIMESTEP) ) {
            diode_real_fft[i] = diode_real[i];
        } else {
            diode_real_fft[i] = 0.;
        }
    }

    // forward FFT
    Tools::realft(diode_real_fft, 1, pad_len * 2);

    // clear previous data as we need new diode response array for new pad_len
    fdiode_real_dynamic_databin.clear(); 

    // save f domain diode response in fdiode_real_dynamic_databin
    for (int i = 0; i < pad_len * 2; i++) {
        fdiode_real_dynamic_databin.push_back( diode_real_fft[i] );
    }
}

In the Setting class, I added the option that control dynamic NFOUR option
Array in the header file
int DYNAMIC_NFOUR;
in the .cc file

DYNAMIC_NFOUR=0;      // 0 : (default) follow pre/user-configured NFOUR value, 1: change NFOUR value by 
following signal length. MK added -2023-05-20-
...
              else if (label == "DYNAMIC_NFOUR") {
                  DYNAMIC_NFOUR = atof( line.substr(line.find_first_of("=") + 1).c_str() );
              }

In the Trigger class, I removed several uncommented lines and the myconvlv_half() and CheckChannelsPass() functions. It looks like we are not using it. But let me know if it has purposes.
In the header file, I removed below arrays

int CheckChannelsPass( vector <double> &V_total_diode);
void myconvlv(vector <double> &data, const int NFOUR, vector <double> &fdiode, vector <double> &diodeconv);
void myconvlv(double *data, const int NFOUR, vector <double> &fdiode, vector <double> &diodeconv);
void myconvlv_half(vector <double> &data, const int NFOUR, vector <double> &fdiode, vector <double> &diodeconv);
void myconvlv_half(double *data, const int NFOUR, vector <double> &fdiode, vector <double> &diodeconv);

in the .cc file, I removed myconvlv_half() and CheckChannelsPass() functions and not used array

iminbin = NFOUR/4 - detector->ibinshift + detector->idelaybeforepeak;
imaxbin = NFOUR/4 - detector->ibinshift + detector->idelaybeforepeak + detector->window;

two myconvlv_half() functions

// test for myconvlv with NFOUR/2 version
void Trigger::myconvlv_half(vector <double> &data,const int NFOUR,vector <double> &fdiode, vector <double> &diodeconv) {
    const int length=NFOUR;
//    double data_copy[length];
    //double fdiode_real[length];
    double power_noise_copy[length];
/*    
    for (int i=0;i<NFOUR/2;i++) {
	data_copy[i]=data[i];
    }
    for (int i=NFOUR/2;i<length;i++) {
	data_copy[i]=0.;
    }
*/  
    for (int i=0;i<length;i++) {
	power_noise_copy[i]=(data[i]*data[i])/Zr*TIMESTEP;
    }
    Tools::realft(power_noise_copy,1,length);
    double ans_copy[length];
    // change the sign (from numerical recipes 648, 649 page)
    for (int j=1;j<length/2;j++) {
	ans_copy[2*j]=(power_noise_copy[2*j]*fdiode[2*j]-power_noise_copy[2*j+1]*fdiode[2*j+1])/((double)length/2);
	//ans_copy[2*j]=(power_noise_copy[2*j]*fdiode[2*j]+power_noise_copy[2*j+1]*fdiode[2*j+1])/((double)length/2);
	ans_copy[2*j+1]=(power_noise_copy[2*j+1]*fdiode[2*j]+power_noise_copy[2*j]*fdiode[2*j+1])/((double)length/2);
	//ans_copy[2*j+1]=(power_noise_copy[2*j+1]*fdiode[2*j]-power_noise_copy[2*j]*fdiode[2*j+1])/((double)length/2);
    }
    ans_copy[0]=power_noise_copy[0]*fdiode[0]/((double)length/2);
    ans_copy[1]=power_noise_copy[1]*fdiode[1]/((double)length/2);
    Tools::realft(ans_copy,-1,length);
    diodeconv.clear();  // remove previous values in diodeconv
    for (int i=0;i<length;i++) {
	//power_noise[i]=power_noise_copy[i];
	//diodeconv[i]=ans_copy[i];
	diodeconv.push_back( ans_copy[i] );
    }
    int iminsamp,imaxsamp; // find min and max samples such that
    // the diode response is fully overlapping with the noise waveform
    iminsamp=(int)(maxt_diode/TIMESTEP);
    // the noise waveform is NFOUR/2 long
    // then for a time maxt_diode/TIMESTEP at the end of that, the
    // diode response function is only partially overlappying with the
    // waveform in the convolution
    imaxsamp=NFOUR/2;
    //  cout << "iminsamp, imaxsamp are " << iminsamp << " " << imaxsamp << "\n";
    if (imaxsamp<iminsamp) {
	cout << "Noise waveform is not long enough for this diode response.\n";
	exit(1);
    }
    //int ibin=((iminsamp+imaxsamp)-(iminsamp+imaxsamp)%2)/2;
    //cout << "ibin is " << ibin << "\n";
    // return the 50th sample, right in the middle
   // onediodeconvl=diodeconv[ibin];
    //mindiodeconvl=0.;
    //for (int i=iminsamp;i<imaxsamp;i++) {
//	if (diodeconv[i]<mindiodeconvl)
//	    mindiodeconvl=diodeconv[i];
  //  }
}

CheckChannelsPass() function

int Trigger::CheckChannelsPass( vector <double> &V_total_diode ) {

    int ch_pass_bin = 0;

    //for (int ibin=iminbin; ibin<imaxbin; ibin++) {
    for (int ibin=(int)(maxt_diode/TIMESTEP);ibin<DATA_BIN_SIZE;ibin++) {

        if ( V_total_diode[ibin] < powerthreshold * rmsdiode ) {
            ch_pass_bin = ibin;
            ibin = DATA_BIN_SIZE;   // stop loop
        }

	    cout << "Noise waveform is not long enough for this diode response.\n";
	    exit(1);
    }

    return ch_pass_bin; // if there was any iwindow values pass the threshold, return = bin number where passed the trigger, if nothing, return = 0;
}

In the Report class, It looks like a lot. But I just removed all... unnecessary repetitive lines. Unfortunately, lots of identical lines are written in the code because It simply has different if conditions. And definitely, It can be shortened to a few lines. I only cleaned (removed repetitive lines) the time-domain signal generation part, triggering part, and readout window setting part
In the header file, I removed functions that loading not-used arrays in the detector class.

void ApplyFilter_NFOUR(int bin_n, Detector *detector, double &vmmhz);
void ApplyPreamp_NFOUR(int bin_n, Detector *detector, double &vmmhz);
void ApplyFOAM_NFOUR(int bin_n, Detector *detector, double &vmmhz);

And I added vector <int> New_NFOUR; array for storing dynamic NFOUR values
in the .cc file,
First, I cleaned the signal generation part in time-domain mode (SIMULATION_MODE == 1)
Mainly, I just cleaned the comment and reorganized the codes.

if (settings1->USE_ARA_ICEATTENU == 1 || settings1->USE_ARA_ICEATTENU == 0)
{
atten_factor = 1. / ray_output[0][ray_sol_cnt] *IceAttenFactor *mag * fresnel;  // assume whichray = 0, now vmmhz1m_tmp has all factors except for the detector properties (antenna gain, etc)
 }
 else if (settings1->USE_ARA_ICEATTENU == 2)
{
atten_factor = 1. / ray_output[0][ray_sol_cnt] *mag * fresnel;  //apply freq dependent IceAttenFactor later
 }

is changed to below

double ice_atten_factor = IceAttenFactor ; ///< assume whichray = 0, now vmmhz1m_tmp has all factors except for the detector properties (antenna gain, etc)
if (settings1->USE_ARA_ICEATTENU == 2) ice_atten_factor = 1.; //apply freq dependent IceAttenFactor later
 atten_factor = 1. / ray_output[0][ray_sol_cnt] *ice_atten_factor *mag * fresnel;

// now new NFOUR for zero padding

 // now we have to make NFOUR/2 number of bins with random init time
//
// as a test, make first as it is and zero pad

double V_forfft[stations[i].strings[j].antennas[k].Nnew[ray_sol_cnt]];
double T_forfft[stations[i].strings[j].antennas[k].Nnew[ray_sol_cnt]];

for (int n = 0; n < stations[i].strings[j].antennas[k].Nnew[ray_sol_cnt]; n++)
 {
// make Tarray, Earray located at the center of Nnew array

T_forfft[n] = Tarray[outbin / 2] - (dT_forfft *(double)(stations[i].strings[j].antennas[k].Nnew[ray_sol_cnt] / 2 - n));
if ((n >= stations[i].strings[j].antennas[k].Nnew[ray_sol_cnt] / 2 - outbin / 2) &&
(n < stations[i].strings[j].antennas[k].Nnew[ray_sol_cnt] / 2 + outbin / 2))
 {
V_forfft[n] = Earray[n - (stations[i].strings[j].antennas[k].Nnew[ray_sol_cnt] / 2 - outbin / 2)];
 }
 else
V_forfft[n] = 0.;

is changed to below

//! second, let's make zero pad what has Nnew length and fill with E-field
int pad_len = stations[i].strings[j].antennas[k].Nnew[ray_sol_cnt];
double V_forfft[pad_len]; ///< E-field value
 double T_forfft[pad_len]; ///< time in ns
double time_offset = Tarray[outbin / 2] - dT_forfft * (double)(pad_len / 2); ///< time offset between E-field center and zero-pad center
int front_pad_len = pad_len / 2 - outbin / 2; ///< length of eampty space in front
int back_pad_len = pad_len / 2 + outbin / 2; ///< length of eampty space in back
for (int n = 0; n < pad_len; n++) {
 //! make Tarray, Earray located at the center of Nnew array
T_forfft[n] = time_offset + (double)n * dT_forfft;
if ((n >= front_pad_len) && (n < back_pad_len)) V_forfft[n] = Earray[n - front_pad_len]; ///< fill with E-field value
else V_forfft[n] = 0.; ///< fill with zero
 }

All the effective height (Heff) lines are simply changed below. I organized values that can be changed by the if condition at the front

int ant_number = 0; 
if (settings1->ANTENNA_MODE == 1) ant_number = k; ///< There are dfferences in Top and Bottom VPol. if ant_number is 2, GetGain_1D_OutZero() will apply Top VPol gain 
int ant_type = detector->stations[i].strings[j].antennas[k].type; ///< 0: VPol, 1: HPol
if (settings1->ALL_ANT_V_ON == 1) ant_type = 0; ///< all the antennas are VPol

//! effective height for last frequency bin 
freq_lastbin = dF_Nnew *((double)pad_len / 2. + 0.5); ///< in Hz 0.5 to place the middle of the bin and avoid zero freq
heff_lastbin = GaintoHeight(detector->GetGain_1D_OutZero(freq_lastbin *1.E-6,   // to MHz
                                                                antenna_theta, antenna_phi, ant_type, ant_number),
                                                                freq_lastbin, icemodel->GetN(detector->stations[i].strings[j].antennas[k]));

//! loop over all frequency bin and apply detector response
 for (int n = 0; n < pad_len / 2; n++) {
//! calculate antenna effective height
    heff = GaintoHeight(detector->GetGain_1D_OutZero(freq_tmp *1.E-6,   // to MHz
                                    antenna_theta, antenna_phi, ant_type, ant_number),
                                    freq_tmp, icemodel->GetN(detector->stations[i].strings[j].antennas[k]));
                                    stations[i].strings[j].antennas[k].Heff[ray_sol_cnt].push_back(heff); ///< save in Report branch

The long trigger scan part is changed to below. I organized values that can be changed by the if condition at the front. So there is no more repeated part. Trust me this is better than before.

                          while (trig_j < ch_ID) {
  
                              int string_i = detector->getStringfromArbAntID(i, trig_j);
                              int antenna_i = detector->getAntennafromArbAntID(i, trig_j);
  
                              //! set channel number for GetThres() and GetThresOffset() at here
                              //! we only use different channle number, if we only want to use VPol for trigger
                              int channel_num = detector->GetChannelfromStringAntenna(i, string_i, antenna_i, settings1);
                              if ((settings1->TRIG_ONLY_LOW_CH_ON == 1) && (settings1->DETECTOR != 3)) channel_num = GetChannelNum8_LowAnt(string_i, antenna_i); ///< reset channel numbers so that bottom antennas have ch 1-8
  
                              //! set threshold value at here
                              //! threshold value can be different hased on NOISE_CHANNEL_MODE and TRIG_ONLY_BH_ON options
                              double Thres_val = detector->GetThres(i, channel_num - 1, settings1) * detector->GetThresOffset(i, channel_num - 1, settings1);
                              if (settings1->NOISE_CHANNEL_MODE == 0) {Thres_val *= trigger->rmsdiode;}
                              else if (settings1->NOISE_CHANNEL_MODE == 1) {Thres_val *= trigger->rmsdiode_ch[channel_num - 1];}
                              else if (settings1->NOISE_CHANNEL_MODE == 2) {
                                  if ((settings1->TRIG_ONLY_BH_ON == 1) && (settings1->DETECTOR == 3)) {Thres_val *= trigger->rmsdiode_ch[channel_num - 1];}
                                  else { 
                                      if (channel_num - 1 < 8) {Thres_val *= trigger->rmsdiode_ch[channel_num - 1];}
                                      else {Thres_val *= trigger->rmsdiode_ch[8];}
                                  }                                
                              }
    
                              //! determine when we launch the trigger search at here 
                              //! if trig_start is true, we do trigger search
                              bool trig_start = false;
                              //! check if we want to use BH chs only for trigger analysis
                              if ((settings1->TRIG_ONLY_BH_ON == 1) && (settings1->DETECTOR == 3)) {
                                  //! trig by BH is only for TestBed case
                                  //! check if this channel is BH ch (DAQchan)
                                  if (detector->stations[i].strings[string_i].antennas[antenna_i].DAQchan == 0) {
                                      trig_start = true;
                                  }
                              }
                              //! in this case just use first 8 chs' thres values
                              else if ((settings1->TRIG_ONLY_LOW_CH_ON == 1) && (settings1->DETECTOR != 3)) {
                                  //! non-TestBed case, and only trig by lower 8 channels
                                  if (antenna_i < 2) {
                                      //! only antenna 0, 1 which are bottom 2 antennas
                                      //! set channel_num as new value (antenna 0, 1 are only possible antennas for channel_num 1 - 8)
                                      trig_start = true;
                                  }
                              }
                              //! other cases, use all possible chs for trigger analysis
                              else {
                                  trig_start = true;
                              }
  
                              //! launch trigger scan
                              if (trig_start == true) {
                                  // cout<<"trig_bin : "<<trig_bin<<endl;
                                  trig_bin = 0;
                                  // cout<<"trig_bin : "<<trig_bin<<endl;
                                  int ant_type = detector->stations[i].strings[string_i].antennas[antenna_i].type;
                                  while (trig_bin < trig_window_bin) {
                                      //! with threshold offset by chs
                                      if (trigger->Full_window[trig_j][trig_i + trig_bin] < Thres_val) {
                                          //! if this channel passed the trigger!
                                          // cout<<"trigger passed at bin "<<trig_i+trig_bin<<" ch : "<<trig_j<<endl;
                                          stations[i].strings[string_i].antennas[antenna_i].Trig_Pass = trig_i + trig_bin;
                                          N_pass++;
                                          if (ant_type == 0) N_pass_V++; // Vpol
                                          if (ant_type == 1) N_pass_H++; // Hpol
                                          if (last_trig_bin < trig_i + trig_bin) last_trig_bin = trig_i + trig_bin;   ///< added for fixed V_mimic
                                          trig_bin = trig_window_bin; ///< if confirmed this channel passed the trigger, no need to do rest of bins
                                          Passed_chs.push_back(trig_j);
                                      }
                                      trig_bin++;
                                  }
                              }
                              //! check all triggered channels not just 3
                              trig_j++;   ///< if station not passed the trigger, just go to next channel
                          }   ///< while trig_j < ch_ID

Lastly I also reorganized readout window part (V_mimic). I organized values that can be changed by the if condition at the front.
terminal msg part is changed to below

//! set msg for triggered event at here
int ant_types = detector->stations[i].strings[string_i].antennas[antenna_i].type;
double posnu_d = event->Nu_Interaction[0].posnu.Distance(detector->stations[i].strings[string_i].antennas[antenna_i]);
 int noise_id = stations[i].strings[string_i].antennas[antenna_i].noise_ID[0];
int likely_sol = stations[i].strings[string_i].antennas[antenna_i].Likely_Sol;
//! skip this passed ch as it already has bin info
if (settings1->TRIG_ONLY_LOW_CH_ON == 0) {
    cout << endl << "trigger passed at bin " << trig_pass << "  passed ch : " << ch_loop << " (" << ant_types << "type) Direct dist btw posnu : " << posnu_d << " noiseID : " << noise_id;
} else if (settings1->TRIG_ONLY_LOW_CH_ON == 1) {
    cout << endl << "trigger passed at bin " << trig_pass << "  passed ant: str[" << string_i << "].ant[" << antenna_i << "] (" << ant_types << "type) Direct dist btw posnu : " << posnu_d << " noiseID : " << noise_id;
}
if (likely_sol != -1) {
    cout << " ViewAngle : " << stations[i].strings[string_i].antennas[antenna_i].view_ang[0] *DEGRAD << " LikelyTrigSignal : " << likely_sol;
}

V_mimic part is changed to below

//! now save the voltage waveform to V_mimic
                                //! set initial time bin index and delay bins length at here
                                int time_bin_i = waveformCenter - waveformLength / 2; ///< center minus half length. yeah, this is initial time
                                int testbed_ch_delay = detector->params.TestBed_Ch_delay_bin[ch_loop] - detector->params.TestBed_BH_Mean_delay_bin; ///< testbed delay
                                int manual_delay_bins = detector->stations[i].strings[string_i].antennas[antenna_i].manual_delay_bin; ///< manual delay by user
                                
                                //! new DAQ waveform writing mechanism test
                                //! loop over readout window bin
                                for (int mimicbin = 0; mimicbin < waveformLength; mimicbin++) {

                                    int time_bin = time_bin_i + mimicbin;
                                    int wf_bin = last_trig_bin + time_bin;
                                    if (V_mimic_mode == 1 || V_mimic_mode == 2) { ///< add delay bin
                                        wf_bin -= testbed_ch_delay;
                                        time_bin -= testbed_ch_delay;
                                    }
                                    if (V_mimic_mode == 2) { ///< add manual delay bin
                                        wf_bin -= manual_delay_bins;
                                        time_bin -= manual_delay_bins;
                                    }
                                    double time_mimic_val = (double)time_bin * settings1->TIMESTEP * 1.e9; ///< set time at here by multiplying time width
                                    if (V_mimic_mode == 2) time_mimic_val += detector->params.TestBed_WFtime_offset_ns; ///< add offset time

                                    stations[i].strings[string_i].antennas[antenna_i].V_mimic.push_back(trigger->Full_window_V[ch_loop][wf_bin] * 1.e3);   // save in mV
                                    stations[i].strings[string_i].antennas[antenna_i].time.push_back(wf_bin);
                                    stations[i].strings[string_i].antennas[antenna_i].time_mimic.push_back(time_mimic_val);    // save in ns

                                    if (mimicbin == 0) {
                                        for (int m = 0; m < stations[i].strings[string_i].antennas[antenna_i].ray_sol_cnt; m++) { ///< calculates time of center of each rays signal based on readout window time config
                                            //! time offset between beginning of readout window and center of signal
                                            double signal_center_offset = (double)(stations[i].strings[string_i].antennas[antenna_i].SignalBin[m] - wf_bin) * settings1->TIMESTEP * 1.e9;
                                            //! time of center of signal based on readout window time config
                                            double signal_center_time = signal_center_offset + time_mimic_val;
                                            stations[i].strings[string_i].antennas[antenna_i].SignalBinTime.push_back(signal_center_time);
                                        }
                                    }
                                }

                                //! set global_trig_bin values
                                //! Global passed bin is the center of the window
                                int global_trig_bins = -1 * time_bin_i; ///< if (V_mimic_mode == 0) Global passed bin is the center of the window
                                if (V_mimic_mode == 1 || V_mimic_mode == 2) global_trig_bins += testbed_ch_delay; ///< Global passed bin is the center of the window + delay to each chs from araGeom
                                if (V_mimic_mode == 2) global_trig_bins += manual_delay_bins; ///< Global passed bin is the center of the window + delay to each chs from araGeom + fitted by eye
                                stations[i].strings[string_i].antennas[antenna_i].global_trig_bin = global_trig_bins;
                            }

Please let me know if you have a comment!