ROI_refinement.cxx

Go to the documentation of this file.

#include "ROI_refinement.h"
#include "PeakFinding.h"
#include <iostream>
#include <set>

using namespace WireCell;
using namespace WireCell::SigProc;

ROI_refinement::ROI_refinement(Waveform::ChannelMaskMap& cmm,int nwire_u, int nwire_v, int nwire_w, float th_factor, float fake_signal_low_th, float fake_signal_high_th, float fake_signal_low_th_ind_factor, float fake_signal_high_th_ind_factor, int pad, int break_roi_loop, float th_peak, float sep_peak, float low_peak_sep_threshold_pre, int max_npeaks, float sigma, float th_percent)
  : nwire_u(nwire_u)
  , nwire_v(nwire_v)
  , nwire_w(nwire_w)
  , th_factor(th_factor)
  , fake_signal_low_th(fake_signal_low_th)
  , fake_signal_high_th(fake_signal_high_th)
  , fake_signal_low_th_ind_factor(fake_signal_low_th_ind_factor)
  , fake_signal_high_th_ind_factor(fake_signal_high_th_ind_factor)
  , pad(pad)
  , break_roi_loop(break_roi_loop)
  , th_peak(th_peak)
  , sep_peak(sep_peak)
  , low_peak_sep_threshold_pre(low_peak_sep_threshold_pre)
  , max_npeaks(max_npeaks)
  , sigma(sigma)
  , th_percent(th_percent)
  , log(Log::logger("sigproc"))
{
  rois_u_tight.resize(nwire_u);
  rois_u_loose.resize(nwire_u);

  rois_v_tight.resize(nwire_v);
  rois_v_loose.resize(nwire_v);

  rois_w_tight.resize(nwire_w);
  
  for (int i=0;i!=nwire_u;i++){
    SignalROIList temp_rois;
    rois_u_tight.at(i) = temp_rois;
  }
  for (int i=0;i!=nwire_v;i++){
    SignalROIList temp_rois;
    rois_v_tight.at(i) = temp_rois;
  }
  for (int i=0;i!=nwire_w;i++){
    SignalROIList temp_rois;
    rois_w_tight.at(i) = temp_rois;
  }
  
  for (int i=0;i!=nwire_u;i++){
    SignalROIList temp_rois;
    rois_u_loose.at(i) = temp_rois;
  }
  for (int i=0;i!=nwire_v;i++){
    SignalROIList temp_rois;
    rois_v_loose.at(i) = temp_rois;
  }

  for (auto it = cmm["bad"].begin(); it!=cmm["bad"].end(); it++){
    int ch = it->first;
    std::vector<std::pair<int,int>> temps;
    bad_ch_map[ch] = temps;
    for (size_t ind = 0; ind < it->second.size(); ind++){
      bad_ch_map[ch].push_back(std::make_pair(it->second.at(ind).first, it->second.at(ind).second));
      //std::cout << ch << " " <<  << std::endl;
    }
  }
  
}

ROI_refinement::~ROI_refinement(){
  Clear();
}

void ROI_refinement::Clear(){
  for (int i=0;i!=nwire_u;i++){
    for (auto it = rois_u_tight.at(i).begin(); it!=rois_u_tight.at(i).end();it++){
      delete *it;
    }
    rois_u_tight.at(i).clear();
    for (auto it = rois_u_loose.at(i).begin(); it!=rois_u_loose.at(i).end();it++){
      delete *it;
    }
    rois_u_loose.at(i).clear();
  }

  for (int i=0;i!=nwire_v;i++){
    for (auto it = rois_v_tight.at(i).begin(); it!=rois_v_tight.at(i).end();it++){
      delete *it;
    }
    rois_v_tight.at(i).clear();
    for (auto it = rois_v_loose.at(i).begin(); it!=rois_v_loose.at(i).end();it++){
      delete *it;
    }
    rois_v_loose.at(i).clear();
  }
  
  for (int i=0;i!=nwire_w;i++){
    for (auto it = rois_w_tight.at(i).begin(); it!=rois_w_tight.at(i).end();it++){
      delete *it;
    }
    rois_w_tight.at(i).clear();
  }
  
  front_rois.clear();
  back_rois.clear();
  contained_rois.clear();
}

void ROI_refinement::apply_roi(int plane, Array::array_xxf& r_data){
  if (plane==0){
    for (int irow = 0 ; irow != r_data.rows(); irow++){
      //refresh ... 
      Waveform::realseq_t signal(r_data.cols(),0);
      // loop ROIs and assign data
      for (auto it = rois_u_loose.at(irow).begin(); it!= rois_u_loose.at(irow).end(); it++){
        SignalROI *roi =  *it;
        // int start_bin = roi->get_start_bin();
        // int end_bin = roi->get_end_bin();

        int start_bin = roi->get_ext_start_bin();
        int end_bin = roi->get_ext_end_bin();
        
        float start_content = r_data(irow,start_bin);
        float end_content = r_data(irow,end_bin);
        for (int i=start_bin; i<end_bin+1; i++){
          int content = r_data(irow,i) - ((end_content - start_content)*(i-start_bin)/(end_bin-start_bin) + start_content);
          signal.at(i) = content;
          //      htemp_signal->SetBinContent(i+1,content);
        }
      }
      // load data back ..
      for (int icol = 0; icol!=r_data.cols();icol++){
        r_data(irow,icol) = signal.at(icol);
      }
    }
  }else if (plane==1){
    for (int irow = 0 ; irow != r_data.rows(); irow++){
      //refresh ... 
      Waveform::realseq_t signal(r_data.cols(),0);
      // loop ROIs and assign data
      for (auto it = rois_v_loose.at(irow).begin(); it!= rois_v_loose.at(irow).end(); it++){
        SignalROI *roi =  *it;
        //int start_bin = roi->get_start_bin();
        //int end_bin = roi->get_end_bin();

        int start_bin = roi->get_ext_start_bin();
        int end_bin = roi->get_ext_end_bin();
        
        float start_content = r_data(irow,start_bin);
        float end_content = r_data(irow,end_bin);
        for (int i=start_bin; i<end_bin+1; i++){
          int content = r_data(irow,i) - ((end_content - start_content)*(i-start_bin)/(end_bin-start_bin) + start_content);
          signal.at(i) = content;
          //htemp_signal->SetBinContent(i+1,content);
        }
      }
      // load data back ..
      for (int icol = 0; icol!=r_data.cols();icol++){
        r_data(irow,icol) = signal.at(icol);
      }
    }
  }else{
    for (int irow = 0 ; irow != r_data.rows(); irow++){
      //refresh ... 
      Waveform::realseq_t signal(r_data.cols(),0);

      // if (irow==69){
      //        std::cout << "Xin: " << rois_w_tight.at(irow).size() << " " << std::endl;
      //        for (auto it = rois_w_tight.at(irow).begin(); it!= rois_w_tight.at(irow).end(); it++){
      //          SignalROI *roi =  *it;
      //          int start_bin = roi->get_ext_start_bin();
      //          int end_bin = roi->get_ext_end_bin();
      //          std::cout << "Xin: " << start_bin << " " << end_bin << std::endl;
      //        }
      // }
      
      // loop ROIs and assign data
      for (auto it = rois_w_tight.at(irow).begin(); it!= rois_w_tight.at(irow).end(); it++){
        SignalROI *roi =  *it;
        //int start_bin = roi->get_start_bin();
        //int end_bin = roi->get_end_bin();

        int start_bin = roi->get_ext_start_bin();
        int end_bin = roi->get_ext_end_bin();
        
        //      float start_content = r_data(irow,start_bin);
        //      float end_content = r_data(irow,end_bin);
        for (int i=start_bin; i<end_bin+1; i++){
          int content = r_data(irow,i);// - ((end_content - start_content)*(i-start_bin)/(end_bin-start_bin) + start_content);
          signal.at(i) = content;
        }
      }
      // load data back ..
      for (int icol = 0; icol!=r_data.cols();icol++){
        r_data(irow,icol) = signal.at(icol);
      }
    }
  }
}

// helper function for getting ROIs
SignalROIChList& ROI_refinement::get_rois_by_plane(int plane)
{
  if (plane==0) return get_u_rois();
  else if (plane==1) return get_v_rois();
  else return get_w_rois();
}

void ROI_refinement::unlink(SignalROI* prev_roi, SignalROI* next_roi){
  if (front_rois.find(prev_roi)!=front_rois.end()){
    SignalROISelection& temp_rois = front_rois[prev_roi];
    auto it = find(temp_rois.begin(),temp_rois.end(),next_roi);
    if (it != temp_rois.end())
      temp_rois.erase(it);
  }
  if (back_rois.find(next_roi)!=back_rois.end()){
    SignalROISelection& temp_rois = back_rois[next_roi];
    auto it = find(temp_rois.begin(),temp_rois.end(),prev_roi);
    if (it != temp_rois.end())
      temp_rois.erase(it);
  }
}

void ROI_refinement::link(SignalROI* prev_roi, SignalROI* next_roi){
  if (front_rois.find(prev_roi)!=front_rois.end()){
    SignalROISelection& temp_rois = front_rois[prev_roi];
    auto it = find(temp_rois.begin(),temp_rois.end(),next_roi);
    if (it == temp_rois.end())
      temp_rois.push_back(next_roi);
  }else{
    SignalROISelection temp_rois;
    temp_rois.push_back(next_roi);
    front_rois[prev_roi] = temp_rois;
  }

  if (back_rois.find(next_roi)!=back_rois.end()){
    SignalROISelection& temp_rois = back_rois[next_roi];
    auto it = find(temp_rois.begin(),temp_rois.end(),prev_roi);
    if (it == temp_rois.end())
      temp_rois.push_back(prev_roi);
  }else{
    SignalROISelection temp_rois;
    temp_rois.push_back(prev_roi);
    back_rois[next_roi] = temp_rois;
  }
}

void ROI_refinement::load_data(int plane, const Array::array_xxf& r_data, ROI_formation& roi_form){
  // fill RMS 
  std::vector<float> plane_rms;
  int offset = 0;
  if (plane==0){
    offset = 0;
    plane_rms = roi_form.get_uplane_rms();
  }else if (plane==1){
    offset = nwire_u;
    plane_rms = roi_form.get_vplane_rms();
  }else if (plane==2){
    offset = nwire_u+nwire_v;
    plane_rms = roi_form.get_wplane_rms();
  }

  // load data ... 
  for (int irow = 0; irow!=r_data.rows(); irow++){
    Waveform::realseq_t signal(r_data.cols());
    if (bad_ch_map.find(irow+offset)!=bad_ch_map.end()){
      for (int icol=0;icol!=r_data.cols();icol++){
        bool flag = true;
        for (size_t i=0; i!=bad_ch_map[irow+offset].size(); i++){
          if (icol >= bad_ch_map[irow+offset].at(i).first &&
              icol <= bad_ch_map[irow+offset].at(i).second){
            flag = false;
            break;
          }
        }
        if (flag){
          signal.at(icol) = r_data(irow,icol);
        }else{
          signal.at(icol) = 0;
        }
      }
    }else{
      for (int icol = 0; icol!= r_data.cols(); icol++){
        signal.at(icol) = r_data(irow,icol);
      }
    }

    int chid = irow+offset;
    // load tight rois
    std::vector<std::pair<int,int>>& uboone_rois = roi_form.get_self_rois(irow+offset);
    for (size_t i=0;i!=uboone_rois.size();i++){
      SignalROI *tight_roi = new SignalROI(plane,irow+offset, uboone_rois.at(i).first,uboone_rois.at(i).second, signal);
      float threshold = plane_rms.at(irow) * th_factor;
      if (tight_roi->get_above_threshold(threshold).size()==0) {
        delete tight_roi;
        continue;
      }
      
      if (plane==0){
        rois_u_tight[chid].push_back(tight_roi);
        if (chid>0){
          //form connectivity map
          for (auto it = rois_u_tight[chid-1].begin();it!=rois_u_tight[chid-1].end();it++){
            SignalROI *prev_roi = *it;
            if (tight_roi->overlap(prev_roi)){
              if (front_rois.find(prev_roi) == front_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(tight_roi);
                front_rois[prev_roi] = temp_rois;
              }else{
                front_rois[prev_roi].push_back(tight_roi);
              }
              if (back_rois.find(tight_roi) == back_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(prev_roi);
                back_rois[tight_roi] = temp_rois;
              }else{
                back_rois[tight_roi].push_back(prev_roi);
              }
            }
          }
        }
        if (chid < nwire_u-1){
          // add the code for the next one to be completed
          for (auto it = rois_u_tight[chid+1].begin();it!=rois_u_tight[chid+1].end();it++){
            SignalROI *next_roi = *it;
            if (tight_roi->overlap(next_roi)){
              if (back_rois.find(next_roi) == back_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(tight_roi);
                back_rois[next_roi] = temp_rois;
              }else{
                back_rois[next_roi].push_back(tight_roi);
              }
              
              if (front_rois.find(tight_roi) == front_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(next_roi);
                front_rois[tight_roi] = temp_rois;
              }else{
                front_rois[tight_roi].push_back(next_roi);
              }
            }
          }
        }
      }else if (plane==1){
        rois_v_tight[chid-nwire_u].push_back(tight_roi);
        if (chid>nwire_u){
          //form connectivity map
          for (auto it = rois_v_tight[chid - nwire_u-1].begin();it!=rois_v_tight[chid-nwire_u-1].end();it++){
            SignalROI *prev_roi = *it;
            if (tight_roi->overlap(prev_roi)){
              if (front_rois.find(prev_roi) == front_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(tight_roi);
                front_rois[prev_roi] = temp_rois;
              }else{
                front_rois[prev_roi].push_back(tight_roi);
              }
              if (back_rois.find(tight_roi) == back_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(prev_roi);
                back_rois[tight_roi] = temp_rois;
              }else{
                back_rois[tight_roi].push_back(prev_roi);
              }
            }
          }
        }
        
        if (chid<nwire_u+nwire_v-1){
          //form connectivity map
          for (auto it = rois_v_tight[chid - nwire_u+1].begin();it!=rois_v_tight[chid-nwire_u+1].end();it++){
            SignalROI *next_roi = *it;
            if (tight_roi->overlap(next_roi)){
              if (back_rois.find(next_roi) == back_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(tight_roi);
                back_rois[next_roi] = temp_rois;
              }else{
                back_rois[next_roi].push_back(tight_roi);
              }
              if (front_rois.find(tight_roi) == front_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(next_roi);
                front_rois[tight_roi] = temp_rois;
              }else{
                front_rois[tight_roi].push_back(next_roi);
              }
            }
          }
        }
      }else {
        rois_w_tight[chid-nwire_u-nwire_v].push_back(tight_roi);
        
        if (chid>nwire_u+nwire_v){
          //form connectivity map
          for (auto it = rois_w_tight[chid-nwire_u-nwire_v-1].begin();it!=rois_w_tight[chid-nwire_u-nwire_v-1].end();it++){
            SignalROI *prev_roi = *it;
            if (tight_roi->overlap(prev_roi)){
              if (front_rois.find(prev_roi) == front_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(tight_roi);
                front_rois[prev_roi] = temp_rois;
              }else{
                front_rois[prev_roi].push_back(tight_roi);
              }
              if (back_rois.find(tight_roi) == back_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(prev_roi);
                back_rois[tight_roi] = temp_rois;
              }else{
                back_rois[tight_roi].push_back(prev_roi);
              }
            }
          }
        }


        if (chid<nwire_u+nwire_v+nwire_w-1){
          //form connectivity map
          for (auto it = rois_w_tight[chid-nwire_u-nwire_v+1].begin();it!=rois_w_tight[chid-nwire_u-nwire_v+1].end();it++){
            SignalROI *next_roi = *it;
            if (tight_roi->overlap(next_roi)){
              if (back_rois.find(next_roi) == back_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(tight_roi);
                back_rois[next_roi] = temp_rois;
              }else{
                back_rois[next_roi].push_back(tight_roi);
              }
              if (front_rois.find(tight_roi) == front_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(next_roi);
                front_rois[tight_roi] = temp_rois;
              }else{
                front_rois[tight_roi].push_back(next_roi);
              }
            }
          }
        }
      }

      
      
    }// loop over tight rois ... 

    if (plane!=2){
      uboone_rois = roi_form.get_loose_rois(chid);
      for (size_t i = 0; i!=uboone_rois.size();i++){
        SignalROI *loose_roi = new SignalROI(plane,chid,uboone_rois.at(i).first,uboone_rois.at(i).second,signal);
        float threshold = plane_rms.at(irow) * th_factor;
        if (loose_roi->get_above_threshold(threshold).size()==0) {
          delete loose_roi;
          continue;
        }
        if (plane==0){
          rois_u_loose[chid].push_back(loose_roi);
          
          if (chid>0){
            //form connectivity map
            for (auto it=rois_u_loose[chid-1].begin();it!=rois_u_loose[chid-1].end();it++){
              SignalROI *prev_roi = *it;
              if (loose_roi->overlap(prev_roi)){
                if (front_rois.find(prev_roi) == front_rois.end()){
                  SignalROISelection temp_rois;
                  temp_rois.push_back(loose_roi);
                  front_rois[prev_roi] = temp_rois;
                }else{
                  front_rois[prev_roi].push_back(loose_roi);
                }
                if (back_rois.find(loose_roi) == back_rois.end()){
                  SignalROISelection temp_rois;
                  temp_rois.push_back(prev_roi);
                  back_rois[loose_roi] = temp_rois;
                }else{
                  back_rois[loose_roi].push_back(prev_roi);
                }
              }
            }
          }
          
          if (chid<nwire_u-1){
            //form connectivity map
            for (auto it=rois_u_loose[chid+1].begin();it!=rois_u_loose[chid+1].end();it++){
              SignalROI *next_roi = *it;
              if (loose_roi->overlap(next_roi)){
                if (back_rois.find(next_roi) == back_rois.end()){
                  SignalROISelection temp_rois;
                  temp_rois.push_back(loose_roi);
                  back_rois[next_roi] = temp_rois;
                }else{
                  back_rois[next_roi].push_back(loose_roi);
                }
                if (front_rois.find(loose_roi) == front_rois.end()){
                  SignalROISelection temp_rois;
                  temp_rois.push_back(next_roi);
                  front_rois[loose_roi] = temp_rois;
                }else{
                  front_rois[loose_roi].push_back(next_roi);
                }
              }
            }
          }
          
          //form contained map ... 
          for (auto it=rois_u_tight[chid].begin();it!=rois_u_tight[chid].end();it++){
            SignalROI *tight_roi = *it;
            if (tight_roi->overlap(loose_roi)){
              if (contained_rois.find(loose_roi)==contained_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(tight_roi);
                contained_rois[loose_roi] = temp_rois;
              }else{
                contained_rois[loose_roi].push_back(tight_roi);
              }
              // if (tight_roi->get_start_bin() < loose_roi->get_start_bin())
              //        std::cout << tight_roi->get_start_bin() << " " << loose_roi->get_start_bin() << " " << tight_roi->get_chid() << " " << tight_roi->get_plane() << std::endl;
            }


          }
        }else if (plane==1){
          rois_v_loose[chid-nwire_u].push_back(loose_roi);
          
          if (chid>nwire_u){
            //form connectivity map
            for (auto it = rois_v_loose[chid-nwire_u-1].begin();it!=rois_v_loose[chid-nwire_u-1].end();it++){
              SignalROI *prev_roi = *it;
              if (loose_roi->overlap(prev_roi)){
                if (front_rois.find(prev_roi) == front_rois.end()){
                  SignalROISelection temp_rois;
                  temp_rois.push_back(loose_roi);
                  front_rois[prev_roi] = temp_rois;
                }else{
                  front_rois[prev_roi].push_back(loose_roi);
                }
                if (back_rois.find(loose_roi) == back_rois.end()){
                  SignalROISelection temp_rois;
                  temp_rois.push_back(prev_roi);
                  back_rois[loose_roi] = temp_rois;
                }else{
                  back_rois[loose_roi].push_back(prev_roi);
                }
              }
            }
          }
          
          if (chid<nwire_u+nwire_v-1){
            //form connectivity map
            for (auto it = rois_v_loose[chid-nwire_u+1].begin();it!=rois_v_loose[chid-nwire_u+1].end();it++){
              SignalROI *next_roi = *it;
              if (loose_roi->overlap(next_roi)){
                if (back_rois.find(next_roi) == back_rois.end()){
                  SignalROISelection temp_rois;
                  temp_rois.push_back(loose_roi);
                  back_rois[next_roi] = temp_rois;
                }else{
                  back_rois[next_roi].push_back(loose_roi);
                }
                if (front_rois.find(loose_roi) == front_rois.end()){
                  SignalROISelection temp_rois;
                  temp_rois.push_back(next_roi);
                  front_rois[loose_roi] = temp_rois;
                }else{
                  front_rois[loose_roi].push_back(next_roi);
                }
              }
            }
          }
          
          //form contained map ... 
          for (auto it = rois_v_tight[chid-nwire_u].begin();it!=rois_v_tight[chid-nwire_u].end();it++){
            SignalROI *tight_roi = *it;
            if (tight_roi->overlap(loose_roi)){
              if (contained_rois.find(loose_roi)==contained_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(tight_roi);
                contained_rois[loose_roi] = temp_rois;
              }else{
                contained_rois[loose_roi].push_back(tight_roi);
              }
            }
          }
        }
      }
    }
    
  } // loop over signal rows
  
}

void ROI_refinement::CleanUpROIs(int plane){

  // clean up ROIs
  std::map<SignalROI*, int> ROIsaved_map;

  if (plane==0){
    //int counter = 0;
    for (size_t i=0;i!=rois_u_loose.size();i++){
      // counter += rois_u_loose.at(i).size();
      for (auto it = rois_u_loose.at(i).begin(); it!= rois_u_loose.at(i).end();it++){
        SignalROI *roi = *it;
        if (ROIsaved_map.find(roi)==ROIsaved_map.end()){
          if (contained_rois.find(roi) != contained_rois.end()){
            // contain good stuff
            SignalROISelection temp_rois;
            temp_rois.push_back(roi);
            ROIsaved_map[roi] = 1;
            
            while(temp_rois.size()){
              SignalROI *temp_roi = temp_rois.back();
              temp_rois.pop_back();
              // save all its neighbour into a temporary holder
              if (front_rois.find(temp_roi)!=front_rois.end()){
                for (auto it1 = front_rois[temp_roi].begin();it1!=front_rois[temp_roi].end();it1++){
                  if (ROIsaved_map.find(*it1)==ROIsaved_map.end()){
                    temp_rois.push_back(*it1);
                    ROIsaved_map[*it1] = 1;
                  }
                }
              }
              if (back_rois.find(temp_roi)!=back_rois.end()){
                for (auto it1 = back_rois[temp_roi].begin();it1!=back_rois[temp_roi].end();it1++){
                  if (ROIsaved_map.find(*it1)==ROIsaved_map.end()){
                    temp_rois.push_back(*it1);
                    ROIsaved_map[*it1] = 1;
                  }
                }
              }
            }
          }
        }
      }
    }
    
    //remove the bad ones ...
    //int counter2 = 0;
    for (size_t i=0;i!=rois_u_loose.size();i++){
      SignalROISelection to_be_removed;
      for (auto it = rois_u_loose.at(i).begin(); it!= rois_u_loose.at(i).end();it++){
        SignalROI *roi = *it;
        if (ROIsaved_map.find(roi) == ROIsaved_map.end()){
          //    counter2 ++;
          to_be_removed.push_back(roi);
          //it = rois_u_loose.at(i).erase(it);
          // check contained map
          if (contained_rois.find(roi)!= contained_rois.end()){
            std::cout << "Wrong! " << std::endl;
          }
          // check front map
          if (front_rois.find(roi)!=front_rois.end()){
            for (auto it1 = front_rois[roi].begin(); it1 != front_rois[roi].end(); it1++){
              auto it2 = find(back_rois[*it1].begin(),back_rois[*it1].end(),roi);
              back_rois[*it1].erase(it2);
            }
            front_rois.erase(roi);
          }
          // check back map
          if (back_rois.find(roi)!=back_rois.end()){
            for (auto it1 = back_rois[roi].begin(); it1!=back_rois[roi].end(); it1++){
              auto it2 = find(front_rois[*it1].begin(),front_rois[*it1].end(),roi);
              front_rois[*it1].erase(it2);
            }
            back_rois.erase(roi);
          }
        }
      }
      
      for (auto it = to_be_removed.begin(); it!= to_be_removed.end(); it++){
        auto it1 = find(rois_u_loose.at(i).begin(), rois_u_loose.at(i).end(),*it);
        rois_u_loose.at(i).erase(it1);
        delete (*it);
      }
    }
  }else if (plane==1){
    
    // int counter = 0;
    for (size_t i=0;i!=rois_v_loose.size();i++){
      //  counter += rois_v_loose.at(i).size();
      for (auto it = rois_v_loose.at(i).begin(); it!= rois_v_loose.at(i).end();it++){
        SignalROI *roi = *it;
        if (ROIsaved_map.find(roi)==ROIsaved_map.end()){
          if (contained_rois.find(roi) != contained_rois.end()){
            // contain good stuff
            SignalROISelection temp_rois;
            temp_rois.push_back(roi);
            ROIsaved_map[roi] = 1;
            
            while(temp_rois.size()){
              SignalROI *temp_roi = temp_rois.back();
              temp_rois.pop_back();
              // save all its neighbour into a temporary holder
              if (front_rois.find(temp_roi)!=front_rois.end()){
                for (auto it1 = front_rois[temp_roi].begin();it1!=front_rois[temp_roi].end();it1++){
                  if (ROIsaved_map.find(*it1)==ROIsaved_map.end()){
                    temp_rois.push_back(*it1);
                    ROIsaved_map[*it1] = 1;
                  }
                }
              }
              if (back_rois.find(temp_roi)!=back_rois.end()){
                for (auto it1 = back_rois[temp_roi].begin();it1!=back_rois[temp_roi].end();it1++){
                  if (ROIsaved_map.find(*it1)==ROIsaved_map.end()){
                    temp_rois.push_back(*it1);
                    ROIsaved_map[*it1] = 1;
                  }
                }
              }
            }
          }
        }
      }
    }
    
    //remove the bad ones ...
    //int counter2 = 0;
    for (size_t i=0;i!=rois_v_loose.size();i++){
      SignalROISelection to_be_removed;
      for (auto it = rois_v_loose.at(i).begin(); it!= rois_v_loose.at(i).end();it++){
        SignalROI *roi = *it;
        if (ROIsaved_map.find(roi) == ROIsaved_map.end()){
          //    counter2 ++;
          to_be_removed.push_back(roi);
          //it = rois_v_loose.at(i).erase(it);
          // check contained map
          if (contained_rois.find(roi)!= contained_rois.end()){
            std::cout << "Wrong! " << std::endl;
          }
          // check front map
          if (front_rois.find(roi)!=front_rois.end()){
            for (auto it1 = front_rois[roi].begin(); it1 != front_rois[roi].end(); it1++){
              auto it2 = find(back_rois[*it1].begin(),back_rois[*it1].end(),roi);
              back_rois[*it1].erase(it2);
            }
            front_rois.erase(roi);
          }
          // check back map
          if (back_rois.find(roi)!=back_rois.end()){
            for (auto it1 = back_rois[roi].begin(); it1!=back_rois[roi].end(); it1++){
              auto it2 = find(front_rois[*it1].begin(),front_rois[*it1].end(),roi);
              front_rois[*it1].erase(it2);
            }
            back_rois.erase(roi);
          }
        }
      }
      
      for (auto it = to_be_removed.begin(); it!= to_be_removed.end(); it++){
        auto it1 = find(rois_v_loose.at(i).begin(), rois_v_loose.at(i).end(),*it);
        rois_v_loose.at(i).erase(it1);
        delete (*it);
      }
    }
  }



  // int counter1 = 0;
  // for (int i=0;i!=rois_v_loose.size();i++){
  //   counter1+=rois_v_loose.at(i).size();
  // }
  
  // std::cout << counter << " " << ROIsaved_map.size() << " " << counter1 << " " << counter2 << std::endl;
}

void ROI_refinement::generate_merge_ROIs(int plane){
  // find tight ROIs not contained by the loose ROIs
  if (plane==0){
    for (int i = 0;i!=nwire_u;i++){
      std::map<SignalROI*,int> covered_tight_rois;
      for (auto it = rois_u_loose.at(i).begin();it!=rois_u_loose.at(i).end();it++){
        SignalROI *roi = *it;
        if (contained_rois.find(roi) != contained_rois.end()){
          for (auto it1 = contained_rois[roi].begin(); it1!= contained_rois[roi].end(); it1++){
            if (covered_tight_rois.find(*it1)==covered_tight_rois.end()){
              covered_tight_rois[*it1]  =1;
            }
          }
        }
      }
      SignalROISelection saved_rois;
      for (auto it = rois_u_tight.at(i).begin();it!=rois_u_tight.at(i).end();it++){
        SignalROI *roi = *it;
        if (covered_tight_rois.find(roi) == covered_tight_rois.end()){
          saved_rois.push_back(roi);
        }
      }
      // if (i == 1212)
      //   std::cout << saved_rois.size() << " " << saved_rois.at(0)->get_start_bin() << " " << saved_rois.at(0)->get_end_bin() << std::endl;
      
      for (auto it = saved_rois.begin(); it!=saved_rois.end();it++){
        SignalROI *roi = *it;
        // Duplicate them 
        SignalROI *loose_roi = new SignalROI(roi);
        
        rois_u_loose.at(i).push_back(loose_roi);
        
        // update all the maps     
        // contained
        SignalROISelection temp_rois;
        temp_rois.push_back(roi);
        contained_rois[loose_roi] = temp_rois;
        // front map loose ROI
        if (i < nwire_u-1){
          for (auto it1 = rois_u_loose.at(i+1).begin(); it1!=rois_u_loose.at(i+1).end(); it1++){
            SignalROI *next_roi = *it1;
            
            if (loose_roi->overlap(next_roi)){
              if (back_rois.find(next_roi) == back_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(loose_roi);
                back_rois[next_roi] = temp_rois;
              }else{
                back_rois[next_roi].push_back(loose_roi);
              }
              
              if (front_rois.find(loose_roi) == front_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(next_roi);
                front_rois[loose_roi] = temp_rois;
              }else{
                front_rois[loose_roi].push_back(next_roi);
              }
            }
            
          }
        }
        // back map loose ROI
        if (i > 0){
          for (auto it1 = rois_u_loose.at(i-1).begin(); it1!=rois_u_loose.at(i-1).end(); it1++){
            SignalROI *prev_roi = *it1;
            if (loose_roi->overlap(prev_roi)){
              if (front_rois.find(prev_roi) == front_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(loose_roi);
                front_rois[prev_roi] = temp_rois;
              }else{
                front_rois[prev_roi].push_back(loose_roi);
              }
              if (back_rois.find(loose_roi) == back_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(prev_roi);
                back_rois[loose_roi] = temp_rois;
              }else{
                back_rois[loose_roi].push_back(prev_roi);
              }
            }
          }
        }
      }
    }
  }else if (plane==1){
    for (int i = 0;i!=nwire_v;i++){
      std::map<SignalROI*,int> covered_tight_rois;
      for (auto it = rois_v_loose.at(i).begin();it!=rois_v_loose.at(i).end();it++){
        SignalROI *roi = *it;
        if (contained_rois.find(roi) != contained_rois.end()){
          for (auto it1 = contained_rois[roi].begin(); it1!= contained_rois[roi].end(); it1++){
            if (covered_tight_rois.find(*it1)==covered_tight_rois.end()){
              covered_tight_rois[*it1]  =1;
            }
          }
        }
      }
      SignalROISelection saved_rois;
      for (auto it = rois_v_tight.at(i).begin();it!=rois_v_tight.at(i).end();it++){
        SignalROI *roi = *it;
        if (covered_tight_rois.find(roi) == covered_tight_rois.end()){
          saved_rois.push_back(roi);
        }
      }
      //if (i == 3885-2400)
      //  std::cout << saved_rois.size() << std::endl;
      //   std::cout << saved_rois.size() << " " << saved_rois.at(0)->get_start_bin() << " " << saved_rois.at(0)->get_end_bin() << std::endl;
      
      for (auto it = saved_rois.begin(); it!=saved_rois.end();it++){
        SignalROI *roi = *it;
        // Duplicate them 
        SignalROI *loose_roi = new SignalROI(roi);
        
        rois_v_loose.at(i).push_back(loose_roi);
        
        // update all the maps     
        // contained
        SignalROISelection temp_rois;
        temp_rois.push_back(roi);
        contained_rois[loose_roi] = temp_rois;
        // front map loose ROI
        if (i < nwire_v-1){
          for (auto it1 = rois_v_loose.at(i+1).begin(); it1!=rois_v_loose.at(i+1).end(); it1++){
            SignalROI *next_roi = *it1;
            
            if (loose_roi->overlap(next_roi)){
              if (back_rois.find(next_roi) == back_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(loose_roi);
                back_rois[next_roi] = temp_rois;
              }else{
                back_rois[next_roi].push_back(loose_roi);
              }
              
              if (front_rois.find(loose_roi) == front_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(next_roi);
                front_rois[loose_roi] = temp_rois;
              }else{
                front_rois[loose_roi].push_back(next_roi);
              }
            }
            
          }
        }
        // back map loose ROI
        if (i > 0){
          for (auto it1 = rois_v_loose.at(i-1).begin(); it1!=rois_v_loose.at(i-1).end(); it1++){
            SignalROI *prev_roi = *it1;
            if (loose_roi->overlap(prev_roi)){
              if (front_rois.find(prev_roi) == front_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(loose_roi);
                front_rois[prev_roi] = temp_rois;
              }else{
                front_rois[prev_roi].push_back(loose_roi);
              }
              if (back_rois.find(loose_roi) == back_rois.end()){
                SignalROISelection temp_rois;
                temp_rois.push_back(prev_roi);
                back_rois[loose_roi] = temp_rois;
              }else{
                back_rois[loose_roi].push_back(prev_roi);
              }
            }
          }
        }
      }
    }
  }
}

void ROI_refinement::CheckROIs(int plane,ROI_formation& roi_form){

  if (plane==0){
    std::vector<float>& rms_u = roi_form.get_uplane_rms();
    
    // for (int i=0;i!=rois_u_loose.size();i++){
    //   for (auto it = rois_u_loose.at(i).begin(); it!= rois_u_loose.at(i).end();it++){
    //     SignalROI *roi = *it;
    //     int chid = roi->get_chid();
    //     if (chid != i) 
    //  std::cout << roi << std::endl;
    
    //     if (front_rois.find(roi)!=front_rois.end()){
    //          for (auto it1 = front_rois[roi].begin();it1!=front_rois[roi].end();it1++){
    //    SignalROI *roi1 = *it1;
    //    int chid1 = roi1->get_chid();
    //    if (chid1!=i+1)
    //      std::cout << roi1 << " " << chid << " " << chid1 << std::endl;
    //  }
    //     }
    //   }
    // }
    
    //loop over u loose
    for (size_t i=0;i!=rois_u_loose.size();i++){
      for (auto it = rois_u_loose.at(i).begin(); it!= rois_u_loose.at(i).end();it++){
        SignalROI *roi = *it;
        int chid = roi->get_chid();
        float th;
        th = th_factor*rms_u.at(chid);
        
        if (front_rois.find(roi)!=front_rois.end()){
          SignalROISelection temp_rois;
          for (auto it1 = front_rois[roi].begin();it1!=front_rois[roi].end();it1++){
            SignalROI *roi1 = *it1;
            int chid1 = roi1->get_chid();
            //std::cout << "F " << i << " " << rois_u_loose.size() << " " << roi1 << " " << chid << " " << chid1 << std::endl;
            float th1;
            th1 = th_factor*rms_u.at(chid1);
            if (roi->overlap(roi1,th,th1)){
            }else{
              temp_rois.push_back(roi1);
              //unlink(roi,roi1);
            }
          }
          for (auto it2 = temp_rois.begin(); it2!= temp_rois.end();it2++){
            unlink(roi,*it2);
          }
        }

        if (back_rois.find(roi)!=back_rois.end()){
          SignalROISelection temp_rois;
          for (auto it1 = back_rois[roi].begin();it1!=back_rois[roi].end();it1++){
            SignalROI *roi1 = *it1;
            int chid1 = roi1->get_chid();
            //std::cout << "B " << roi1 << " " << chid << " " << chid1 << std::endl;
            float th1;
            th1 = th_factor*rms_u.at(chid1);
            if (roi->overlap(roi1,th,th1)){
            }else{
              temp_rois.push_back(roi1);
              //unlink(roi,roi1);
            }
          }
          for (auto it2 = temp_rois.begin(); it2!= temp_rois.end();it2++){
            unlink(*it2,roi);
          }
        }
        
      }
    }
  }else if (plane==1){
    std::vector<float>& rms_v = roi_form.get_vplane_rms();
    //loop over v loose
    for (size_t i=0;i!=rois_v_loose.size();i++){
      for (auto it = rois_v_loose.at(i).begin(); it!= rois_v_loose.at(i).end();it++){
        SignalROI *roi = *it;
        int chid = roi->get_chid()-nwire_u;
        float th;
        th = th_factor*rms_v.at(chid);
        if (front_rois.find(roi)!=front_rois.end()){
          SignalROISelection temp_rois;
          for (auto it1 = front_rois[roi].begin();it1!=front_rois[roi].end();it1++){
            SignalROI *roi1 = *it1;
            int chid1 = roi1->get_chid()-nwire_u;
            float th1;
            th1 = th_factor*rms_v.at(chid1);
            if (roi->overlap(roi1,th,th1)){
            }else{
              temp_rois.push_back(roi1);
              // unlink(roi,roi1);
            }
          }
          for (auto it2 = temp_rois.begin(); it2!= temp_rois.end();it2++){
            unlink(roi,*it2);
          }
        }
        
        if (back_rois.find(roi)!=back_rois.end()){
          SignalROISelection temp_rois;
          for (auto it1 = back_rois[roi].begin();it1!=back_rois[roi].end();it1++){
            SignalROI *roi1 = *it1;
            int chid1 = roi1->get_chid()-nwire_u;
            float th1;
            th1 = th_factor*rms_v.at(chid1);
            if (roi->overlap(roi1,th,th1)){
            }else{
              temp_rois.push_back(roi1);
              // unlink(roi,roi1);
            }
          }
          for (auto it2 = temp_rois.begin(); it2!= temp_rois.end();it2++){
            unlink(*it2,roi);
          }
        }
      }
    }
  }
}

void ROI_refinement::CleanUpCollectionROIs(){
  // deal with tight ROIs, 
  // scan with all the tight ROIs to look for peaks above certain threshold, put in a temporary set
  float mean_threshold = fake_signal_low_th;
  float threshold = fake_signal_high_th; //electrons, about 1/2 of MIP per tick ...
  std::set<SignalROI*> Good_ROIs;
  for (int i=0;i!=nwire_w;i++){
    for (auto it = rois_w_tight.at(i).begin();it!=rois_w_tight.at(i).end();it++){
      SignalROI* roi = *it;
      //std::cout << "Xin: " << roi->get_max_height() << " " << roi->get_average_heights() << std::endl;
      if (roi->get_above_threshold(threshold).size()!=0 || roi->get_average_heights() > mean_threshold)
        Good_ROIs.insert(roi);
    }
  }
  // for a particular ROI if it is not in, or it is not connected with one in the temporary map, then remove it
  std::list<SignalROI*> Bad_ROIs;
  for (int i=0;i!=nwire_w;i++){
    for (auto it = rois_w_tight.at(i).begin();it!=rois_w_tight.at(i).end();it++){
      SignalROI* roi = *it;
      
      if (Good_ROIs.find(roi)!=Good_ROIs.end()) continue;
      if (front_rois.find(roi)!=front_rois.end()){
        SignalROISelection next_rois = front_rois[roi];
        int flag_qx = 0;
        for (size_t i=0;i!=next_rois.size();i++){
          SignalROI* roi1 = next_rois.at(i);
          if (Good_ROIs.find(roi1)!=Good_ROIs.end()) {
            flag_qx = 1;
            continue;
          }
        }
        if (flag_qx == 1) continue;
      }
      
      if (back_rois.find(roi)!=back_rois.end()){
        SignalROISelection next_rois = back_rois[roi];
        int flag_qx = 0;
        for (size_t i=0;i!=next_rois.size();i++){
          SignalROI* roi1 = next_rois.at(i);
          if (Good_ROIs.find(roi1)!=Good_ROIs.end()) {
            flag_qx = 1;
            continue;
          }
        }
        if (flag_qx == 1) continue;
      }
      
      Bad_ROIs.push_back(roi);
    }
  }
  
  // remove the ROI and then update the map
  
  for (auto it = Bad_ROIs.begin(); it!=Bad_ROIs.end(); it ++){
    SignalROI* roi = *it;
    int chid = roi->get_chid()-nwire_u-nwire_v;
    //std::cout << chid << std::endl;
    if (front_rois.find(roi)!=front_rois.end()){
      SignalROISelection next_rois = front_rois[roi];
      for (size_t i=0;i!=next_rois.size();i++){
        //unlink the current roi
        unlink(roi,next_rois.at(i));
      }
      front_rois.erase(roi);
    }
    
    if (back_rois.find(roi)!=back_rois.end()){
      SignalROISelection prev_rois = back_rois[roi];
      for (size_t i=0;i!=prev_rois.size();i++){
        //unlink the current roi
        unlink(prev_rois.at(i),roi);
      }
      back_rois.erase(roi);
    }
    auto it1 = find(rois_w_tight.at(chid).begin(), rois_w_tight.at(chid).end(),roi);
    if (it1 != rois_w_tight.at(chid).end())
      rois_w_tight.at(chid).erase(it1);
    
    delete roi;
  }
  
}

void ROI_refinement::CleanUpInductionROIs(int plane){
   // deal with loose ROIs
  // focus on the isolated ones first
  float mean_threshold = fake_signal_low_th;
  float threshold = fake_signal_high_th;
  mean_threshold *= fake_signal_low_th_ind_factor;
  threshold *= fake_signal_high_th_ind_factor;
  
  std::list<SignalROI*> Bad_ROIs;
  if (plane==0){
    for (int i=0;i!=nwire_u;i++){
      for (auto it = rois_u_loose.at(i).begin();it!=rois_u_loose.at(i).end();it++){
        SignalROI* roi = *it;
        //std::cout << "Xin: " << roi->get_max_height() << " " << roi->get_average_heights() << std::endl;
        if (front_rois.find(roi)==front_rois.end() && back_rois.find(roi)==back_rois.end()){
          if (roi->get_above_threshold(threshold).size()==0 && roi->get_average_heights() < mean_threshold)
            Bad_ROIs.push_back(roi);
        }
      }
    }
    for (auto it = Bad_ROIs.begin(); it!=Bad_ROIs.end(); it ++){
      SignalROI* roi = *it;
      int chid = roi->get_chid();
      auto it1 = find(rois_u_loose.at(chid).begin(), rois_u_loose.at(chid).end(),roi);
      if (it1 != rois_u_loose.at(chid).end())
        rois_u_loose.at(chid).erase(it1);

      if (front_rois.find(roi)!=front_rois.end()){
        SignalROISelection next_rois = front_rois[roi];
        for (size_t i=0;i!=next_rois.size();i++){
          //unlink the current roi
          unlink(roi,next_rois.at(i));
        }
        front_rois.erase(roi);
      }
      
      if (back_rois.find(roi)!=back_rois.end()){
        SignalROISelection prev_rois = back_rois[roi];
        for (size_t i=0;i!=prev_rois.size();i++){
          //unlink the current roi
          unlink(prev_rois.at(i),roi);
        }
        back_rois.erase(roi);
      }
      
      delete roi;
    }
    Bad_ROIs.clear();
  }else if (plane==1){
    for (int i=0;i!=nwire_v;i++){
      for (auto it = rois_v_loose.at(i).begin();it!=rois_v_loose.at(i).end();it++){
        SignalROI* roi = *it;
        //      std::cout << "Xin: " << roi->get_max_height() << " " << roi->get_average_heights() << std::endl;
        if (front_rois.find(roi)==front_rois.end() && back_rois.find(roi)==back_rois.end()){
          if (roi->get_above_threshold(threshold).size()==0 && roi->get_average_heights() < mean_threshold)
            Bad_ROIs.push_back(roi);
        }
      }
    }
    for (auto it = Bad_ROIs.begin(); it!=Bad_ROIs.end(); it ++){
      SignalROI* roi = *it;
      int chid = roi->get_chid()-nwire_u;
      auto it1 = find(rois_v_loose.at(chid).begin(), rois_v_loose.at(chid).end(),roi);
      if (it1 != rois_v_loose.at(chid).end())
        rois_v_loose.at(chid).erase(it1);

      if (front_rois.find(roi)!=front_rois.end()){
        SignalROISelection next_rois = front_rois[roi];
        for (size_t i=0;i!=next_rois.size();i++){
          //unlink the current roi
          unlink(roi,next_rois.at(i));
        }
        front_rois.erase(roi);
      }
      
      if (back_rois.find(roi)!=back_rois.end()){
        SignalROISelection prev_rois = back_rois[roi];
        for (size_t i=0;i!=prev_rois.size();i++){
          //unlink the current roi
          unlink(prev_rois.at(i),roi);
        }
        back_rois.erase(roi);
      }
      
      
      delete roi;
    }
  }


  // threshold = fake_signal_low_th;
  std::set<SignalROI*> Good_ROIs;
  if (plane==0){
    for (int i=0;i!=nwire_u;i++){
      for (auto it = rois_u_loose.at(i).begin();it!=rois_u_loose.at(i).end();it++){
        SignalROI* roi = *it;
        if (roi->get_above_threshold(threshold).size()!=0 || roi->get_average_heights() > mean_threshold)
          Good_ROIs.insert(roi);
      }
    }
    Bad_ROIs.clear();
    for (int i=0;i!=nwire_u;i++){
      for (auto it = rois_u_loose.at(i).begin();it!=rois_u_loose.at(i).end();it++){
        SignalROI* roi = *it;
        
        if (Good_ROIs.find(roi)!=Good_ROIs.end()) continue;
        if (front_rois.find(roi)!=front_rois.end()){
          SignalROISelection next_rois = front_rois[roi];
          int flag_qx = 0;
          for (size_t i=0;i!=next_rois.size();i++){
            SignalROI* roi1 = next_rois.at(i);
            if (Good_ROIs.find(roi1)!=Good_ROIs.end()) {
              flag_qx = 1;
              continue;
            }
          }
          if (flag_qx == 1) continue;
        }
        
        if (back_rois.find(roi)!=back_rois.end()){
          SignalROISelection next_rois = back_rois[roi];
          int flag_qx = 0;
          for (size_t i=0;i!=next_rois.size();i++){
            SignalROI* roi1 = next_rois.at(i);
            if (Good_ROIs.find(roi1)!=Good_ROIs.end()) {
              flag_qx = 1;
              continue;
            }
          }
          if (flag_qx == 1) continue;
        }
        
        Bad_ROIs.push_back(roi);
      }
    }
    for (auto it = Bad_ROIs.begin(); it!=Bad_ROIs.end(); it ++){
      SignalROI* roi = *it;
      int chid = roi->get_chid();
      //std::cout << chid << std::endl;
      if (front_rois.find(roi)!=front_rois.end()){
        SignalROISelection next_rois = front_rois[roi];
        for (size_t i=0;i!=next_rois.size();i++){
          //unlink the current roi
          unlink(roi,next_rois.at(i));
        }
        front_rois.erase(roi);
      }
      
      if (back_rois.find(roi)!=back_rois.end()){
        SignalROISelection prev_rois = back_rois[roi];
        for (size_t i=0;i!=prev_rois.size();i++){
          //unlink the current roi
          unlink(prev_rois.at(i),roi);
        }
        back_rois.erase(roi);
      }
      auto it1 = find(rois_u_loose.at(chid).begin(), rois_u_loose.at(chid).end(),roi);
      if (it1 != rois_u_loose.at(chid).end())
        rois_u_loose.at(chid).erase(it1);
      
      delete roi;
    }
  }else if (plane==1){
    
    Good_ROIs.clear();
    for (int i=0;i!=nwire_v;i++){
      for (auto it = rois_v_loose.at(i).begin();it!=rois_v_loose.at(i).end();it++){
        SignalROI* roi = *it;
        if (roi->get_above_threshold(threshold).size()!=0)
          Good_ROIs.insert(roi);
      }
    }
    Bad_ROIs.clear();
    for (int i=0;i!=nwire_v;i++){
      for (auto it = rois_v_loose.at(i).begin();it!=rois_v_loose.at(i).end();it++){
        SignalROI* roi = *it;
        
        if (Good_ROIs.find(roi)!=Good_ROIs.end()) continue;
        if (front_rois.find(roi)!=front_rois.end()){
          SignalROISelection next_rois = front_rois[roi];
          int flag_qx = 0;
          for (size_t i=0;i!=next_rois.size();i++){
            SignalROI* roi1 = next_rois.at(i);
            if (Good_ROIs.find(roi1)!=Good_ROIs.end()) {
              flag_qx = 1;
              continue;
            }
          }
          if (flag_qx == 1) continue;
        }
        
        if (back_rois.find(roi)!=back_rois.end()){
          SignalROISelection next_rois = back_rois[roi];
          int flag_qx = 0;
          for (size_t i=0;i!=next_rois.size();i++){
            SignalROI* roi1 = next_rois.at(i);
            if (Good_ROIs.find(roi1)!=Good_ROIs.end()) {
              flag_qx = 1;
              continue;
            }
          }
          if (flag_qx == 1) continue;
        }
        
        Bad_ROIs.push_back(roi);
      }
    }
    for (auto it = Bad_ROIs.begin(); it!=Bad_ROIs.end(); it ++){
      SignalROI* roi = *it;
      int chid = roi->get_chid()-nwire_u;
      //std::cout << chid << std::endl;
      if (front_rois.find(roi)!=front_rois.end()){
        SignalROISelection next_rois = front_rois[roi];
        for (size_t i=0;i!=next_rois.size();i++){
          //unlink the current roi
          unlink(roi,next_rois.at(i));
        }
        front_rois.erase(roi);
      }
      
      if (back_rois.find(roi)!=back_rois.end()){
        SignalROISelection prev_rois = back_rois[roi];
        for (size_t i=0;i!=prev_rois.size();i++){
          //unlink the current roi
          unlink(prev_rois.at(i),roi);
        }
        back_rois.erase(roi);
      }
      auto it1 = find(rois_v_loose.at(chid).begin(), rois_v_loose.at(chid).end(),roi);
      if (it1 != rois_v_loose.at(chid).end())
        rois_v_loose.at(chid).erase(it1);
      
      delete roi;
    }
  }
}

void ROI_refinement::ShrinkROI(SignalROI *roi, ROI_formation& roi_form){
  // get tight ROI as a inner boundary
  // get the nearby ROIs with threshold as some sort of boundary 
  int start_bin = roi->get_start_bin();
  int end_bin = roi->get_end_bin();
  if (start_bin <0 || end_bin <0) return;
  
  int chid = roi->get_chid();
  int plane = roi->get_plane();
  std::vector<float>& contents = roi->get_contents();
  
  float threshold1=0;
  if (plane==0){
    threshold1 = roi_form.get_uplane_rms().at(chid) * th_factor;
  }else if (plane==1){
    threshold1 = roi_form.get_vplane_rms().at(chid-nwire_u)* th_factor;
  }
  
  int channel_save = 1240;
  int print_flag = 0;

  // std::cout << "check tight ROIs " << std::endl;
  // use to save contents
  Waveform::realseq_t temp_signal(end_bin-start_bin+1,0);
  // TH1F *htemp = new TH1F("htemp","htemp",end_bin-start_bin+1,start_bin,end_bin+1);
  
  // check tight ROIs
  if (contained_rois.find(roi)!=contained_rois.end()){
    for (auto it = contained_rois[roi].begin();it!=contained_rois[roi].end();it++){
      SignalROI *tight_roi = *it;
      int start_bin1 = tight_roi->get_start_bin();
      int end_bin1 = tight_roi->get_end_bin();
      
      if (chid == channel_save && print_flag)
        std::cout << "Tight "  " " << start_bin1 << " " << end_bin1 << std::endl;

      for (int i=start_bin1;i<=end_bin1;i++){
        if (i-start_bin >=0 && i-start_bin <int(temp_signal.size())){
          //      htemp->SetBinContent(i-start_bin+1,1);
          temp_signal.at(i-start_bin) = 1;
        }
      }
    }
  }

  // std::cout << "check front ROIs " << std::endl;

  //check front ROIs
  if (front_rois.find(roi)!=front_rois.end()){
    for (auto it=front_rois[roi].begin();it!=front_rois[roi].end();it++){
      SignalROI *next_roi = *it;
      int start_bin1 = next_roi->get_start_bin();
      int chid1 = next_roi->get_chid();
      int plane1 = next_roi->get_plane();
      float threshold=0;
      if (plane1==0){
        threshold = roi_form.get_uplane_rms().at(chid1) * th_factor;
      }else if (plane1==1){
        threshold = roi_form.get_vplane_rms().at(chid1-nwire_u) * th_factor;
      }
      std::vector<std::pair<int,int>> contents_above_threshold = next_roi->get_above_threshold(threshold);
      for (size_t i=0;i!=contents_above_threshold.size();i++){
        if (chid == channel_save && print_flag)
          std::cout << "Front " << chid1 << " " << start_bin1 + contents_above_threshold.at(i).first << " " << start_bin1 + contents_above_threshold.at(i).second << std::endl;

        for (int j=contents_above_threshold.at(i).first;j<=contents_above_threshold.at(i).second;j++){
          if (j+start_bin1-start_bin >=0 && j+start_bin1-start_bin < int(temp_signal.size())){
            if (contents.at(j+start_bin1-start_bin) > threshold1)
              temp_signal.at(j+start_bin1-start_bin) = 1;
            //        htemp->SetBinContent(j+start_bin1-start_bin+1,1);
          }
        }
      }
    }
  }
  
  //std::cout << "check back ROIs " << std::endl;

  //check back ROIs
  if (back_rois.find(roi)!=back_rois.end()){
    for (auto it=back_rois[roi].begin();it!=back_rois[roi].end();it++){
      SignalROI *prev_roi = *it;
      int start_bin1 = prev_roi->get_start_bin();
      int chid1 = prev_roi->get_chid();
      int plane1 = prev_roi->get_plane();
      float threshold = 0;
      if (plane1==0){
        threshold = roi_form.get_uplane_rms().at(chid1) * th_factor;
      }else if (plane1==1){
        threshold = roi_form.get_vplane_rms().at(chid1-nwire_u)* th_factor;
      }
      std::vector<std::pair<int,int>> contents_above_threshold = prev_roi->get_above_threshold(threshold);
      for (size_t i=0;i!=contents_above_threshold.size();i++){
        if (chid == channel_save && print_flag)
          std::cout << "Back " << chid1 << " " << start_bin1 + contents_above_threshold.at(i).first << " " << start_bin1 + contents_above_threshold.at(i).second << std::endl;

        for (int j=contents_above_threshold.at(i).first;j<=contents_above_threshold.at(i).second;j++){
          if (j+start_bin1-start_bin >=0 && j+start_bin1-start_bin <int(temp_signal.size())){
            if (contents.at(j+start_bin1-start_bin) > threshold1)
              temp_signal.at(j+start_bin1-start_bin) = 1;
            //        htemp->SetBinContent(j+start_bin1-start_bin+1,1);
          }
        }
      }
    }
  }

  //std::cout << "check contents " << std::endl;

  // // consider the 1/2 of the peak as threshold;
  // float max = 0;
  // for (int i=0;i!=contents.size();i++){
  //   if (contents.at(i) > max)
  //     max = contents.at(i);
  // }
  // for (int i=0;i!=contents.size();i++){
  //   // if (contents.at(i) > max/2. && contents.at(i) > threshold1*2 ) htemp->SetBinContent(i+1,1);
  // }
  
  // get the first bin, and last bin, add pad
  int new_start_bin=start_bin;
  int new_end_bin=end_bin;
  for (size_t i=0;i!= temp_signal.size(); i++){
    if (temp_signal.at(i) >0){
      new_start_bin = i + start_bin;
      break;
    }
  }
  for (int i = int(temp_signal.size())-1;i>=0;i--){
    if (temp_signal.at(i) > 0){
      new_end_bin = i + start_bin;
      break;
    }
  }
  new_start_bin -= pad;
  new_end_bin += pad;
  if (new_start_bin < start_bin) new_start_bin = start_bin;
  if (new_end_bin > end_bin) new_end_bin = end_bin;
  
  if (chid == channel_save && print_flag)
    std::cout << "check contents " << " " << start_bin << " " << end_bin << " " << new_start_bin << " " << new_end_bin << std::endl;
  

  // create a new ROI
  Waveform::realseq_t signal(end_bin+1);
  // TH1F *h1 = new TH1F("h1","h1",end_bin+1,0,end_bin+1);
  for (int i=new_start_bin; i<=new_end_bin;i++){
    signal.at(i) = contents.at(i-start_bin);
    //   h1->SetBinContent(i+1,contents.at(i-start_bin));
  }
  
  SignalROISelection new_rois;
  if (new_start_bin >=0 && new_end_bin > new_start_bin){
    SignalROI *new_roi = new SignalROI(plane,chid,new_start_bin,new_end_bin,signal);
    new_rois.push_back(new_roi);
  }

  // std::cout << "update maps " << std::endl;
  
  // update the list 
  if (chid < nwire_u){
    auto it = std::find(rois_u_loose.at(chid).begin(),rois_u_loose.at(chid).end(),roi);
    rois_u_loose.at(chid).erase(it);
    for (size_t i=0;i!=new_rois.size();i++){
      rois_u_loose.at(chid).push_back(new_rois.at(i));
    }
  }else if (chid < nwire_u+nwire_v){
    auto it = std::find(rois_v_loose.at(chid-nwire_u).begin(),rois_v_loose.at(chid-nwire_u).end(),roi);
    rois_v_loose.at(chid-nwire_u).erase(it);
    for (size_t i=0;i!=new_rois.size();i++){
      rois_v_loose.at(chid-nwire_u).push_back(new_rois.at(i));
    }
  }
  
  // update all the maps 
  // update front map
  if (front_rois.find(roi)!=front_rois.end()){
    SignalROISelection next_rois = front_rois[roi];
    for (size_t i=0;i!=next_rois.size();i++){
      //unlink the current roi
      unlink(roi,next_rois.at(i));
      //loop new rois and link them
      for (size_t j=0;j!=new_rois.size();j++){
        if (new_rois.at(j)->overlap(next_rois.at(i)))
          link(new_rois.at(j),next_rois.at(i));
      }
    }
    front_rois.erase(roi);
  }
  // update back map
  if (back_rois.find(roi)!=back_rois.end()){
    SignalROISelection prev_rois = back_rois[roi];
    for (size_t i=0;i!=prev_rois.size();i++){
      // unlink the current roi
      unlink(prev_rois.at(i),roi);
      // loop new rois and link them
      for (size_t j=0;j!=new_rois.size();j++){
        if (new_rois.at(j)->overlap(prev_rois.at(i)))
          link(prev_rois.at(i),new_rois.at(j));
      }
    }
    back_rois.erase(roi);
  }
  
  // update contained map 
  if (contained_rois.find(roi)!=contained_rois.end()){
    SignalROISelection tight_rois = contained_rois[roi];
    for (size_t i=0;i!=tight_rois.size();i++){
      for (size_t j=0;j!=new_rois.size();j++){
        if (new_rois.at(j)->overlap(tight_rois.at(i))){
          if (contained_rois.find(new_rois.at(j)) == contained_rois.end()){
            SignalROISelection temp_rois;
            temp_rois.push_back(tight_rois.at(i));
            contained_rois[new_rois.at(j)] = temp_rois;
          }else{
            contained_rois[new_rois.at(j)].push_back(tight_rois.at(i));
          }
        }
      }
    }
    contained_rois.erase(roi);
  }
  
  // delete the old ROI
  delete roi;

  // delete htemp;
  // delete h1;
}

void ROI_refinement::ShrinkROIs(int plane, ROI_formation& roi_form){
  // collect all ROIs
  SignalROISelection all_rois;
  if (plane==0){
    for (size_t i=0;i!=rois_u_loose.size();i++){
      for (auto it = rois_u_loose.at(i).begin(); it!= rois_u_loose.at(i).end(); it++){
        all_rois.push_back(*it);
      }
    }
  }else if (plane==1){
    for (size_t i=0;i!=rois_v_loose.size();i++){
      for (auto it = rois_v_loose.at(i).begin(); it!= rois_v_loose.at(i).end(); it++){
        all_rois.push_back(*it);
      }
    }
  }
  for (size_t i=0;i!=all_rois.size();i++){
    ShrinkROI(all_rois.at(i),roi_form);
  }
}

void ROI_refinement::BreakROI(SignalROI *roi, float rms){

  //  std::cout << "haha " << std::endl;
  
  // main algorithm 
  int start_bin = roi->get_start_bin();
  int end_bin = roi->get_end_bin();
  
  if (start_bin <0 || end_bin <0 ) return;

  // if (roi->get_chid()==1240){
  //   std::cout << "xin: " << start_bin << " " << end_bin << std::endl;
  // }

  Waveform::realseq_t temp_signal(end_bin-start_bin+1,0);
  // TH1F *htemp = new TH1F("htemp","htemp",end_bin-start_bin+1,start_bin,end_bin+1);
  std::vector<float>& contents = roi->get_contents();
  for (size_t i=0;i!=temp_signal.size();i++){
    temp_signal.at(i) = contents.at(i);
    //    htemp->SetBinContent(i+1,contents.at(i));
  }
  
  // float th_peak = 3.0;
  // float sep_peak = 6.0;
  float low_peak_sep_threshold = low_peak_sep_threshold_pre; // electrons
  if (low_peak_sep_threshold < sep_peak * rms) 
    low_peak_sep_threshold = sep_peak * rms;
  std::set<int> saved_boundaries;

  PeakFinding s(max_npeaks, sigma, th_percent);
  const int nfound = s.find_peak(temp_signal);
  // TSpectrum *s = new TSpectrum(200);
  // Int_t nfound = s->Search(htemp,2,"nobackground new",0.1);

  if (nfound == max_npeaks) {
    log->debug("ROI_refinement: local ch index {} (plane {}), found max peaks {} with threshold={}",
               roi->get_chid(), roi->get_plane(), nfound, th_percent);
  }
  
  if (nfound > 1){
    int npeaks = s.GetNPeaks();
    double *peak_pos = s.GetPositionX();
    double  *peak_height = s.GetPositionY();
    //const int temp_length = max_npeaks + 5;
    std::vector<int> order_peak_pos;
    int npeaks_threshold = 0;
    for (int j=0;j!=npeaks;j++){
      order_peak_pos.push_back(*(peak_pos+j) + start_bin);
      if (*(peak_height+j)>th_peak*rms){
        npeaks_threshold ++;
      }
    }

    
    
    if (npeaks_threshold >1){
      std::sort(order_peak_pos.begin(),order_peak_pos.end());
      float valley_pos[205];
      valley_pos[0] = start_bin;

      // find the first real valley
      float min = 1e9;
      for (int k=0; k< order_peak_pos[0]-start_bin;k++){
        if (temp_signal.at(k) < min){
          min = temp_signal.at(k);
          valley_pos[0] = k+start_bin;
        }
      }
      if (valley_pos[0] != start_bin){
        // temporarily push n a value ...
        order_peak_pos.push_back(0);
        //
        for (int j=npeaks-1;j>=0;j--){
          order_peak_pos[j+1] = order_peak_pos[j];
        }
        npeaks ++;
        order_peak_pos[0] = start_bin;
        for (int j=start_bin; j!=valley_pos[0];j++){
          if (temp_signal.at(j-start_bin) > temp_signal.at(order_peak_pos[0]-start_bin))
            //        htemp->GetBinContent(j-start_bin+1) > htemp->GetBinContent(order_peak_pos[0]-start_bin+1))
            order_peak_pos[0] = j;
        }
        valley_pos[0] = start_bin;
      }

      //   std::cout << "kaka1 " << npeaks << std::endl;
      
      for (int j=0;j!=npeaks-1;j++){
        float min = 1e9;
        valley_pos[j+1] = order_peak_pos[j];

        //      std::cout << order_peak_pos[j]-start_bin << " " << order_peak_pos[j+1]-start_bin << std::endl;
        
        for (int k = order_peak_pos[j]-start_bin; k< order_peak_pos[j+1]-start_bin;k++){
          if (temp_signal.at(k) < min){
            min = temp_signal.at(k);
            valley_pos[j+1] = k+start_bin;
          }
        }
      }
      
      // std::cout << "kaka1 " << npeaks << std::endl;
      
      //find the end ... 
      valley_pos[npeaks] = end_bin;
      min = 1e9;
      for (int k=order_peak_pos[npeaks-1]-start_bin; k<= end_bin-start_bin;k++){
        if (temp_signal.at(k) < min){ 
          min = temp_signal.at(k);
          valley_pos[npeaks] = k+start_bin;
        }
      }

      //  std::cout << "kaka1 " << npeaks << std::endl;
      
      if (valley_pos[npeaks]!=end_bin){
        npeaks ++;
        valley_pos[npeaks] = end_bin;
        // temporarily add an elements ...
        order_peak_pos.push_back(end_bin);
        //      order_peak_pos[npeaks-1] = end_bin;
        for (int j=valley_pos[npeaks-1];j!=valley_pos[npeaks];j++){
          if (temp_signal.at(j-start_bin) > temp_signal.at(order_peak_pos[npeaks-1] -start_bin))
            order_peak_pos[npeaks-1] = j;
        }
      }

      // std::cout << "kaka1 " << npeaks << std::endl;
      
      // if (roi->get_chid() == 1195 && roi->get_plane() == WirePlaneType_t(0)){
      //        for (int j=0;j!=npeaks;j++){
      //          std::cout << valley_pos[j] << " " << htemp->GetBinContent(valley_pos[j]-start_bin+1)<< " " << order_peak_pos[j] << " " << htemp->GetBinContent( order_peak_pos[j]-start_bin+1) << " " << valley_pos[j+1] << " " << htemp->GetBinContent(valley_pos[j+1] - start_bin+1)<< " " << rms * sep_peak << std::endl ;
      //        }
      // }
        

      
      // need to organize the peaks and valleys ... 
      float valley_pos1[205];
      float peak_pos1[205];
      int npeaks1 = 0;
      // fill in the first valley;
      valley_pos1[0] = valley_pos[0];
      for (int j=0;j<npeaks;j++){
        if (npeaks1 >0){
          //std::cout << valley_pos[j]-start_bin << " " << valley_pos1[npeaks1]-start_bin << std::endl;
          // find the lowest valley except the first peak, except the first one
          if (temp_signal.at(valley_pos[j]-start_bin) < temp_signal.at(valley_pos1[npeaks1]-start_bin)){
            valley_pos1[npeaks1] = valley_pos[j];
          }
        }
        
        // if (roi->get_chid() == 1195 && roi->get_plane() == WirePlaneType_t(0)){
        //   std::cout << "c: " << order_peak_pos[j] << " " << htemp->GetBinContent(order_peak_pos[j]-start_bin+1) << " " << valley_pos1[npeaks1] << " " << htemp->GetBinContent(valley_pos1[npeaks1]-start_bin+1) << std::endl;
        // }

        // find the next peak
        if (temp_signal.at(order_peak_pos[j]-start_bin) - temp_signal.at(valley_pos1[npeaks1]-start_bin) > low_peak_sep_threshold){
          peak_pos1[npeaks1] = order_peak_pos[j] ;
          npeaks1 ++;
          int flag1 = 0;

          for (int k=j+1;k!=npeaks+1;k++){
            // find the highest peak before ... 
            if (k<=npeaks){
              if (temp_signal.at(order_peak_pos[k-1]-start_bin) > temp_signal.at(peak_pos1[npeaks1-1]-start_bin))
                peak_pos1[npeaks1-1] = order_peak_pos[k-1];
            }
            
            if (temp_signal.at(peak_pos1[npeaks1-1]-start_bin) - temp_signal.at(valley_pos[k]-start_bin) > low_peak_sep_threshold){
              valley_pos1[npeaks1] = valley_pos[k];
              j = k-1;
              flag1 = 1;
              break;
            }
            // find the next valley
          }
          if (flag1 == 0){
            valley_pos1[npeaks1] = valley_pos[npeaks];
            j = npeaks;
          }

          // if (roi->get_chid() == 1240 && roi->get_plane() == 0)
          //   std::cout << "c: " << npeaks << " " << valley_pos1[npeaks1-1] << " " << peak_pos1[npeaks1-1] << " " << valley_pos1[npeaks1] << " " << rms * sep_peak << std::endl;
        }
      }
      // fill the last valley
      valley_pos1[npeaks1] = valley_pos[npeaks];
      //std::cout << roi->get_plane() << " " << roi->get_chid() << " " << npeaks << " " << npeaks1 << " " ;
      //std::cout << start_bin << " " << end_bin << std::endl;


      // if (roi->get_chid() == 1240 && roi->get_plane() == 0){
      //        for (int j=0;j!=npeaks1;j++){
      //          std::cout << valley_pos1[j] << " " << peak_pos1[j] << " " <<  valley_pos1[j+1] << std::endl ;
      //        }
      // }
      

      //      std::cout << "kaka1 " << std::endl;

      for (int j=0;j!=npeaks1;j++){
        int start_pos = valley_pos1[j];
        int end_pos = valley_pos1[j+1];
        
        // if (roi->get_chid()==1195)
        //   std::cout << "b " << start_pos << " " << end_pos << std::endl;
        
        saved_boundaries.insert(start_pos);
        saved_boundaries.insert(end_pos);
      }

      Waveform::realseq_t temp1_signal = temp_signal;
      temp_signal.clear();
      temp_signal.resize(temp1_signal.size(),0);
      //      htemp->Reset();
      for (int j=0;j!=npeaks1;j++){
        //int flag = 0;
        int start_pos = valley_pos1[j];
        double start_content = temp1_signal.at(valley_pos1[j]-start_bin); //htemp1->GetBinContent(valley_pos1[j]-start_bin+1);
        int end_pos = valley_pos1[j+1];
        double end_content = temp1_signal.at(valley_pos1[j+1]-start_bin);//htemp1->GetBinContent(valley_pos1[j+1]-start_bin+1);
        
        //      std::cout << "a " << start_pos << " " << end_pos << std::endl;

        if (saved_boundaries.find(start_pos) != saved_boundaries.end() ||
            saved_boundaries.find(end_pos) != saved_boundaries.end()){
          // if (roi->get_chid() == 1240 && roi->get_plane() == 0)
          //   std::cout << "d: " << start_pos << " " << end_pos << std::endl;

          for (int k = start_pos; k!=end_pos+1;k++){
            double temp_content = temp1_signal.at(k-start_bin) -  (start_content + (end_content-start_content) * (k-start_pos) / (end_pos - start_pos));
            temp_signal.at(k-start_bin) = temp_content;
          }
        }
      }
      //     delete htemp1;

      //std::cout << "kaka2 " << std::endl;
      // loop through the tight ROIs it contains ...
      // if nothing in the range, add things back ... 
      // if (roi->get_chid() == 1151){
      //        std::cout << "Break:  "  << roi->get_chid() << " " << start_bin << " " << end_bin << " " << nfound << std::endl;
      //        std::cout << npeaks1 << std::endl;
      //        std::cout << contained_rois[roi].size() << std::endl;
      
      for (auto it = contained_rois[roi].begin(); it!= contained_rois[roi].end(); it++){
        SignalROI *temp_roi = *it;
        int temp_flag = 0;
        for (int i=temp_roi->get_start_bin(); i<= temp_roi->get_end_bin(); i++){
          //std::cout << i-roi->get_start_bin() << " " << i << " " << temp_roi->get_start_bin() << " " << temp_roi->get_end_bin() << " " << roi->get_start_bin() << " " << roi->get_end_bin() << " " << roi->get_chid() << " " << roi->get_plane() << std::endl;
          if (i-int(roi->get_start_bin())>=0 && i-int(roi->get_start_bin()) < int(temp_signal.size()))
            if (temp_signal.at(i-roi->get_start_bin())!=0){
              temp_flag = 1;
              break;
            }
        }
        //      std::cout << temp_flag << std::endl;
        if (temp_flag==0){
          for (int i=temp_roi->get_start_bin(); i<= temp_roi->get_end_bin(); i++){
            // std::cout << i-roi->get_start_bin() << " " << i-temp_roi->get_start_bin() << std::endl;
             if (i-int(roi->get_start_bin())>=0 && i-int(roi->get_start_bin()) < int(temp_signal.size()))
              temp_signal.at(i-roi->get_start_bin()) = temp_roi->get_contents().at(i-temp_roi->get_start_bin());
            //      htemp->SetBinContent(i-roi->get_start_bin()+1, );
          }
        }
        
      }
      // } // if (1151)

    }
  }

  //std::cout << "kaka3 " << std::endl;
  
  
  // if (roi->get_chid() == 1151){
  //   std::cout << "Break:  "  << roi->get_chid() << " " << start_bin << " " << end_bin << " " << nfound << std::endl;
  //   for (int i=0;i<htemp->GetNbinsX();i++){
  //     std::cout << htemp->GetBinContent(i+1) <<std::endl;
  //   }
  // }

  for (int qx = 0; qx!=2; qx++){
    // Now we should go through the system again and re-adjust the content
    std::vector<std::pair<int,int>> bins;
    for (int i=0;i<int(temp_signal.size());i++){
      if (temp_signal.at(i) < th_factor*rms){
        int start = i;
        int end = i;
        for (int j=i+1;j< int(temp_signal.size());j++){
          if (temp_signal.at(j) < th_factor*rms){
            end = j;
          }else{
            break;
          }
        }
        bins.push_back(std::make_pair(start,end));
        // if (roi->get_chid() == 1240)
        //   std::cout << qx << " " <<  start+start_bin << " " << end + start_bin << " " << 3*rms  << std::endl;
        i = end;
      }
    }

    // std::cout << "kaka4 " << std::endl;

    std::vector<int> saved_b;
    for (int i=0;i!=int(bins.size());i++){
      int start = bins.at(i).first;
      int end = bins.at(i).second;
      // find minimum or zero
      float min = 1e9;
      int bin_min = start;
      for (int j=start;j<=end;j++){
        if (fabs(temp_signal.at(j)) < 1e-3){
          bin_min = j;
          break;
        }else{
          if (temp_signal.at(j) < min){
            min = temp_signal.at(j); 
            bin_min = j;
          }
        }
      }
      
      // if (roi->get_chid() == 1308)
      //        std::cout << bin_min+start_bin << std::endl;
      
      saved_b.push_back(bin_min);
    }

    // std::cout << "kaka5 " << std::endl;
    // test
    
    //    if (saved_b.size() >=0){
    {
      Waveform::realseq_t temp1_signal = temp_signal;
      temp_signal.clear();
      temp_signal.resize(temp1_signal.size(),0);
      //      htemp->Reset();
      // std::cout << saved_b.size() << " " << bins.size() << " " << htemp->GetNbinsX() << std::endl;
      for (int j=0;j<int(saved_b.size())-1;j++){
        // if (irow==1240)
        //   std::cout << saved_b.size() << " " << j << " " << saved<< std::endl;

        
        //int flag = 0;
        int start_pos = saved_b.at(j);
        float start_content = temp1_signal.at(start_pos);//htemp1->GetBinContent(start_pos+1);
        int end_pos = saved_b.at(j+1);
        float end_content = temp1_signal.at(end_pos);//htemp1->GetBinContent(end_pos+1);
        
        for (int k = start_pos; k!=end_pos+1;k++){
          double temp_content = temp1_signal.at(k) - (start_content + (end_content-start_content) * (k-start_pos) / (end_pos - start_pos));
          temp_signal.at(k) = temp_content;
        }
      }
    }
  }

  // std::cout << "kaka6 " << std::endl;
  
  // get back to the original content 
  for (int i=0;i!=int(temp_signal.size());i++){
    contents.at(i) = temp_signal.at(i);//htemp->GetBinContent(i+1);
  }
  
  //  delete s;
  //  delete htemp;
   
}

void ROI_refinement::BreakROI1(SignalROI *roi){
  int start_bin = roi->get_start_bin();
  int end_bin = roi->get_end_bin();
  if (start_bin <0 || end_bin < 0) return;

  Waveform::realseq_t temp_signal(end_bin-start_bin+1,0);
  // TH1F *htemp = new TH1F("htemp","htemp",end_bin-start_bin+1,start_bin,end_bin+1);
  std::vector<float>& contents = roi->get_contents();
  for (size_t i=0;i!=temp_signal.size();i++){
    temp_signal.at(i) = contents.at(i);
    //   htemp->SetBinContent(i+1,contents.at(i));
  }

  // now create many new ROIs
  std::vector<int> bins;
  for (size_t i=0;i!= temp_signal.size(); i++){
    if (fabs(temp_signal.at(i))<1e-3)
      bins.push_back(i+start_bin);
  }
  int chid = roi->get_chid();
  int plane = roi->get_plane();
  SignalROISelection new_rois;

  // if (chid == 1274)
  //   std::cout << "BreakROI1: " << chid << " " << roi->get_start_bin() << " " << roi->get_end_bin() << " " << bins.size()  << " " << htemp->GetBinContent(1) << " " << htemp->GetBinContent(end_bin-start_bin+1) << std::endl;

  Waveform::realseq_t signal(end_bin+1,0);
  //  TH1F *h1 = new TH1F("h1","h1",end_bin+1,0,end_bin+1);
  for (int i=0;i!=int(bins.size())-1;i++){
    int start_bin1 = bins.at(i);
    int end_bin1 = bins.at(i+1);
    // if (chid == 1274)
    //   std::cout << start_bin1 << " " << end_bin1 << std::endl;
    signal.clear();
    signal.resize(end_bin+1,0);
    //    h1->Reset();
    for (int j=start_bin1;j<=end_bin1;j++){
      signal.at(j) = temp_signal.at(j-start_bin);
      //      h1->SetBinContent(j+1,htemp->GetBinContent(j-start_bin+1));
    }
    if (start_bin1 >=0 && end_bin1 >start_bin1){
      SignalROI *sub_roi = new SignalROI(plane,chid,start_bin1,end_bin1,signal);
      new_rois.push_back(sub_roi);
    }
  }

  // update the list 
  if (chid < nwire_u){
    auto it = std::find(rois_u_loose.at(chid).begin(),rois_u_loose.at(chid).end(),roi);
    rois_u_loose.at(chid).erase(it);
    for (size_t i=0;i!=new_rois.size();i++){
      rois_u_loose.at(chid).push_back(new_rois.at(i));
    }
  }else if (chid < nwire_u+nwire_v){
    auto it = std::find(rois_v_loose.at(chid-nwire_u).begin(),rois_v_loose.at(chid-nwire_u).end(),roi);
    rois_v_loose.at(chid-nwire_u).erase(it);
    for (size_t i=0;i!=new_rois.size();i++){
      rois_v_loose.at(chid-nwire_u).push_back(new_rois.at(i));
    }
  }
  
  // update all the maps 
  // update front map
  if (front_rois.find(roi)!=front_rois.end()){
    SignalROISelection next_rois = front_rois[roi];
    for (size_t i=0;i!=next_rois.size();i++){
      //unlink the current roi
      unlink(roi,next_rois.at(i));
      //loop new rois and link them
      for (size_t j=0;j!=new_rois.size();j++){
        if (new_rois.at(j)->overlap(next_rois.at(i)))
          link(new_rois.at(j),next_rois.at(i));
      }
    }
    front_rois.erase(roi);
  }
  // update back map
  if (back_rois.find(roi)!=back_rois.end()){
    SignalROISelection prev_rois = back_rois[roi];
    for (size_t i=0;i!=prev_rois.size();i++){
      // unlink the current roi
      unlink(prev_rois.at(i),roi);
      // loop new rois and link them
      for (size_t j=0;j!=new_rois.size();j++){
        if (new_rois.at(j)->overlap(prev_rois.at(i)))
          link(prev_rois.at(i),new_rois.at(j));
      }
    }
    back_rois.erase(roi);
  }
  
  // update contained map 
  if (contained_rois.find(roi)!=contained_rois.end()){
    SignalROISelection tight_rois = contained_rois[roi];
    for (size_t i=0;i!=tight_rois.size();i++){
      for (size_t j=0;j!=new_rois.size();j++){
        if (new_rois.at(j)->overlap(tight_rois.at(i))){
          if (contained_rois.find(new_rois.at(j)) == contained_rois.end()){
            SignalROISelection temp_rois;
            temp_rois.push_back(tight_rois.at(i));
            contained_rois[new_rois.at(j)] = temp_rois;
          }else{
            contained_rois[new_rois.at(j)].push_back(tight_rois.at(i));
          }
        }
      }
    }
    contained_rois.erase(roi);
  }
  
  // delete the old ROI
  delete roi;
  //  delete h1;
  //  delete htemp;
}

void ROI_refinement::BreakROIs(int plane, ROI_formation& roi_form){
  SignalROISelection all_rois;

  if (plane==0){
    std::vector<float>& rms_u = roi_form.get_uplane_rms();
    for (size_t i=0;i!=rois_u_loose.size();i++){
      for (auto it = rois_u_loose.at(i).begin(); it!= rois_u_loose.at(i).end(); it++){
        
        BreakROI(*it,rms_u.at(i));
        all_rois.push_back(*it);
        
      }
    }
  }else if (plane==1){
    std::vector<float>& rms_v = roi_form.get_vplane_rms();
    for (size_t i=0;i!=rois_v_loose.size();i++){
      for (auto it = rois_v_loose.at(i).begin(); it!= rois_v_loose.at(i).end(); it++){
        BreakROI(*it,rms_v.at(i));
        all_rois.push_back(*it);
      }
    }
  }

  
  for (size_t i=0;i!=all_rois.size();i++){
    // if (all_rois.at(i)->get_chid()==1151){
    //   std::cout << all_rois.at(i)->get_chid() << " " << all_rois.at(i)->get_start_bin() << " " << all_rois.at(i)->get_end_bin() << std::endl;
    //   for (int j=0;j!=all_rois.at(i)->get_contents().size();j++){
    //  std::cout << j << " " << all_rois.at(i)->get_contents().at(j) << std::endl;  
    //   }
    // }
    
    BreakROI1(all_rois.at(i));
  }

  
}


void ROI_refinement::refine_data(int plane, ROI_formation& roi_form){

  //if (plane==2) std::cout << "Xin: " << rois_w_tight.at(69).size() << " " << std::endl;
  
  //std::cout << "Clean up loose ROIs" << std::endl;
  CleanUpROIs(plane);
  //std::cout << "Generate more loose ROIs from isolated good tight ROIs" << std::endl;
  generate_merge_ROIs(plane);

  // if (plane==2)  std::cout << "Xin: " << rois_w_tight.at(69).size() << " " << std::endl;
  
  for (int qx = 0; qx!=break_roi_loop; qx++){
    // std::cout << "Break loose ROIs" << std::endl;
    BreakROIs(plane, roi_form);
    // std::cout << "Clean up ROIs 2nd time" << std::endl;
    CheckROIs(plane, roi_form);
    CleanUpROIs(plane);
  }

  // if (plane==2)  std::cout << "Xin: " << rois_w_tight.at(69).size() << " " << std::endl;
  
  
  //  std::cout << "Shrink ROIs" << std::endl;
  ShrinkROIs(plane, roi_form);
  // std::cout << "Clean up ROIs 3rd time" << std::endl;
  CheckROIs(plane, roi_form);
  CleanUpROIs(plane);

  //  if (plane==2)  std::cout << "Xin: " << rois_w_tight.at(69).size() << " " << std::endl;

  // Further reduce fake hits
  // std::cout << "Remove fake hits " << std::endl;
  if (plane==2){
    CleanUpCollectionROIs();
  }else{
    CleanUpInductionROIs(plane);
  }

  ExtendROIs();
  //TestROIs();
  //if (plane==2)  std::cout << "Xin: " << rois_w_tight.at(69).size() << " " << std::endl;
}

// Basically rewrite the refine_data()
// But now do the ROI operation step by step
void ROI_refinement::refine_data_debug_mode(int plane, ROI_formation& roi_form, const std::string& cmd){

  if (cmd=="CleanUpROIs") {
    // std::cout << "Clean up loose ROIs" << std::endl;
    CleanUpROIs(plane);
    //std::cout << "Generate more loose ROIs from isolated good tight ROIs" << std::endl;
    generate_merge_ROIs(plane);
  }

  else if (cmd=="BreakROIs") {
    // std::cout << "Break loose ROIs" << std::endl;
    BreakROIs(plane, roi_form);
    // std::cout << "Clean up ROIs 2nd time" << std::endl;
    CheckROIs(plane, roi_form);
    CleanUpROIs(plane);
  }

  else if (cmd=="ShrinkROIs") {
    //  std::cout << "Shrink ROIs" << std::endl;
    ShrinkROIs(plane, roi_form);
    // std::cout << "Clean up ROIs 3rd time" << std::endl;
    CheckROIs(plane, roi_form);
    CleanUpROIs(plane);
  }

  else if (cmd=="ExtendROIs") {
    // Further reduce fake hits
    // std::cout << "Remove fake hits " << std::endl;
    if (plane==2){
      CleanUpCollectionROIs();
    }else{
      CleanUpInductionROIs(plane);
    }

    ExtendROIs();
    //TestROIs();
  }

 }

void ROI_refinement::TestROIs(){
  
  for (int chid = 0; chid != nwire_u; chid ++){
    for (auto it = rois_u_loose.at(chid).begin(); it!= rois_u_loose.at(chid).end();it++){
      SignalROI *roi =  *it;
      //loop through front
      for (auto it1 = front_rois[roi].begin(); it1!=front_rois[roi].end(); it1++){
        SignalROI *roi1 = *it1;
        if (find(rois_u_loose.at(chid+1).begin(), rois_u_loose.at(chid+1).end(), roi1) == rois_u_loose.at(chid+1).end())
          std::cout << chid << " u " << +1 << " " << roi << " " << roi1 << std::endl;
      }
        
      // loop through back 
      for (auto it1 = back_rois[roi].begin(); it1!=back_rois[roi].end(); it1++){
        SignalROI *roi1 = *it1;
        if (find(rois_u_loose.at(chid-1).begin(), rois_u_loose.at(chid-1).end(), roi1) == rois_u_loose.at(chid-1).end())
          std::cout << chid << " u " << -1 << " " << roi << " " << roi1 << std::endl;
      }
    }
  }


   for (int chid = 0; chid != nwire_v; chid ++){
    for (auto it = rois_v_loose.at(chid).begin(); it!= rois_v_loose.at(chid).end();it++){
      SignalROI *roi =  *it;
      //loop through front
      for (auto it1 = front_rois[roi].begin(); it1!=front_rois[roi].end(); it1++){
        SignalROI *roi1 = *it1;
        if (find(rois_v_loose.at(chid+1).begin(), rois_v_loose.at(chid+1).end(), roi1) == rois_v_loose.at(chid+1).end())
          std::cout << chid << " v " << +1 << " " << roi << " " << roi1 << std::endl;
      }
        
      // loop through back 
      for (auto it1 = back_rois[roi].begin(); it1!=back_rois[roi].end(); it1++){
        SignalROI *roi1 = *it1;
        if (find(rois_v_loose.at(chid-1).begin(), rois_v_loose.at(chid-1).end(), roi1) == rois_v_loose.at(chid-1).end())
          std::cout << chid << " v " << -1 << " " << roi << " " << roi1 << std::endl;
      }
    }
  }
   
   for (int chid = 0; chid != nwire_w; chid ++){
     for (auto it = rois_w_tight.at(chid).begin(); it!= rois_w_tight.at(chid).end();it++){
       SignalROI *roi =  *it;
       //loop through front
       for (auto it1 = front_rois[roi].begin(); it1!=front_rois[roi].end(); it1++){
         SignalROI *roi1 = *it1;
        if (find(rois_w_tight.at(chid+1).begin(), rois_w_tight.at(chid+1).end(), roi1) == rois_w_tight.at(chid+1).end())
          std::cout << chid << " w " << +1 << " " << roi << " " << roi1 << std::endl;
      }
        
      // loop through back 
      for (auto it1 = back_rois[roi].begin(); it1!=back_rois[roi].end(); it1++){
        SignalROI *roi1 = *it1;
        if (find(rois_w_tight.at(chid-1).begin(), rois_w_tight.at(chid-1).end(), roi1) == rois_w_tight.at(chid-1).end())
          std::cout << chid << " w " << -1 << " " << roi << " " << roi1 << std::endl;
      }
    }
  }
}


void ROI_refinement::ExtendROIs(){

  bool flag = true;
  
  // U plane ... 
  for (int chid = 0; chid != nwire_u; chid ++){
    rois_u_loose.at(chid).sort(CompareRois());
    for (auto it = rois_u_loose.at(chid).begin(); it!= rois_u_loose.at(chid).end();it++){
      SignalROI *roi =  *it;
      // initialize the extended bins ... 
      roi->set_ext_start_bin(roi->get_start_bin());
      roi->set_ext_end_bin(roi->get_end_bin());

      if (flag){
        //loop through front
        for (auto it1 = front_rois[roi].begin(); it1!=front_rois[roi].end(); it1++){
          SignalROI *roi1 = *it1;
          int ext_start_bin = roi->get_ext_start_bin();
          int ext_end_bin = roi->get_ext_end_bin();
          if (ext_start_bin > roi1->get_start_bin()) ext_start_bin = roi1->get_start_bin();
          if (ext_end_bin < roi1->get_end_bin()) ext_end_bin = roi1->get_end_bin();
          roi->set_ext_start_bin(ext_start_bin);
          roi->set_ext_end_bin(ext_end_bin);
          //    std::cout << roi->get_chid() << " " << roi1->get_chid() << " " << roi->get_start_bin() << " " << roi->get_end_bin() << " " << roi1->get_start_bin() << " " << roi1->get_end_bin() << std::endl;
        }
        
        // loop through back 
        for (auto it1 = back_rois[roi].begin(); it1!=back_rois[roi].end(); it1++){
          SignalROI *roi1 = *it1;
          int ext_start_bin = roi->get_ext_start_bin();
          int ext_end_bin = roi->get_ext_end_bin();
          if (ext_start_bin > roi1->get_start_bin()) ext_start_bin = roi1->get_start_bin();
          if (ext_end_bin < roi1->get_end_bin()) ext_end_bin = roi1->get_end_bin();
          roi->set_ext_start_bin(ext_start_bin);
          roi->set_ext_end_bin(ext_end_bin);
          //    std::cout << roi->get_chid() << " " << roi1->get_chid() << " " << roi->get_start_bin() << " " << roi->get_end_bin() << " " << roi1->get_start_bin() << " " << roi1->get_end_bin() << std::endl;
        }
        //std::cout << chid << " " << roi->get_start_bin() << " " << roi->get_end_bin() << std::endl;
      }
    }

    if (flag){
      SignalROI *prev_roi = 0;
      for (auto it = rois_u_loose.at(chid).begin(); it!= rois_u_loose.at(chid).end();it++){
        SignalROI *roi =  *it;
        if (prev_roi!=0){
          if (prev_roi->get_ext_end_bin() > roi->get_ext_start_bin()){
            prev_roi->set_ext_end_bin(int(prev_roi->get_end_bin() * 0.5 + roi->get_start_bin()*0.5));
            roi->set_ext_start_bin(int(prev_roi->get_end_bin() * 0.5 + roi->get_start_bin()*0.5));
          }
        }
        prev_roi = roi;
      }
    }
    // for (auto it = rois_u_loose.at(chid).begin(); it!= rois_u_loose.at(chid).end();it++){
    //   SignalROI *roi =  *it;
    //   std::cout << chid << " " << roi->get_ext_start_bin() << " " << roi->get_ext_end_bin() << std::endl;
    // }
  }
  // V plane
  for (int chid = 0; chid != nwire_v; chid ++){
    rois_v_loose.at(chid).sort(CompareRois());
    for (auto it = rois_v_loose.at(chid).begin(); it!= rois_v_loose.at(chid).end();it++){
      SignalROI *roi =  *it;
      // initialize the extended bins ... 
      roi->set_ext_start_bin(roi->get_start_bin());
      roi->set_ext_end_bin(roi->get_end_bin());

      if (flag){
        //loop through front
        for (auto it1 = front_rois[roi].begin(); it1!=front_rois[roi].end(); it1++){
          SignalROI *roi1 = *it1;
          int ext_start_bin = roi->get_ext_start_bin();
          int ext_end_bin = roi->get_ext_end_bin();
          if (ext_start_bin > roi1->get_start_bin()) ext_start_bin = roi1->get_start_bin();
          if (ext_end_bin < roi1->get_end_bin()) ext_end_bin = roi1->get_end_bin();
          roi->set_ext_start_bin(ext_start_bin);
          roi->set_ext_end_bin(ext_end_bin);
          //    std::cout << roi->get_chid() << " " << roi1->get_chid() << " " << roi->get_start_bin() << " " << roi->get_end_bin() << " " << roi1->get_start_bin() << " " << roi1->get_end_bin() << std::endl;
        }
        
        // loop through back 
        for (auto it1 = back_rois[roi].begin(); it1!=back_rois[roi].end(); it1++){
          SignalROI *roi1 = *it1;
          int ext_start_bin = roi->get_ext_start_bin();
          int ext_end_bin = roi->get_ext_end_bin();
          if (ext_start_bin > roi1->get_start_bin()) ext_start_bin = roi1->get_start_bin();
          if (ext_end_bin < roi1->get_end_bin()) ext_end_bin = roi1->get_end_bin();
          roi->set_ext_start_bin(ext_start_bin);
          roi->set_ext_end_bin(ext_end_bin);
          //    std::cout << roi->get_chid() << " " << roi1->get_chid() << " " << roi->get_start_bin() << " " << roi->get_end_bin() << " " << roi1->get_start_bin() << " " << roi1->get_end_bin() << std::endl;
        }
        //std::cout << chid << " " << roi->get_start_bin() << " " << roi->get_end_bin() << std::endl;
      }
    }

    if (flag){
      SignalROI *prev_roi = 0;
      for (auto it = rois_v_loose.at(chid).begin(); it!= rois_v_loose.at(chid).end();it++){
        SignalROI *roi =  *it;
        if (prev_roi!=0){
          if (prev_roi->get_ext_end_bin() > roi->get_ext_start_bin()){
            prev_roi->set_ext_end_bin(int(prev_roi->get_end_bin() * 0.5 + roi->get_start_bin()*0.5));
            roi->set_ext_start_bin(int(prev_roi->get_end_bin() * 0.5 + roi->get_start_bin()*0.5));
          }
        }
        prev_roi = roi;
      }
    }
    // for (auto it = rois_v_loose.at(chid).begin(); it!= rois_v_loose.at(chid).end();it++){
    //   SignalROI *roi =  *it;
    //   std::cout << chid << " " << roi->get_ext_start_bin() << " " << roi->get_ext_end_bin() << std::endl;
    // }
  }

  // W plane
  for (int chid = 0; chid != nwire_w; chid ++){
    rois_w_tight.at(chid).sort(CompareRois());
    for (auto it = rois_w_tight.at(chid).begin(); it!= rois_w_tight.at(chid).end();it++){
      SignalROI *roi =  *it;
      // initialize the extended bins ... 
      roi->set_ext_start_bin(roi->get_start_bin());
      roi->set_ext_end_bin(roi->get_end_bin());

      if (flag){
        //loop through front
        for (auto it1 = front_rois[roi].begin(); it1!=front_rois[roi].end(); it1++){
          SignalROI *roi1 = *it1;
          int ext_start_bin = roi->get_ext_start_bin();
          int ext_end_bin = roi->get_ext_end_bin();
          if (ext_start_bin > roi1->get_start_bin()) ext_start_bin = roi1->get_start_bin();
          if (ext_end_bin < roi1->get_end_bin()) ext_end_bin = roi1->get_end_bin();
          roi->set_ext_start_bin(ext_start_bin);
          roi->set_ext_end_bin(ext_end_bin);
          //    std::cout << roi->get_chid() << " " << roi1->get_chid() << " " << roi->get_start_bin() << " " << roi->get_end_bin() << " " << roi1->get_start_bin() << " " << roi1->get_end_bin() << std::endl;
        }
        
        // loop through back 
        for (auto it1 = back_rois[roi].begin(); it1!=back_rois[roi].end(); it1++){
          SignalROI *roi1 = *it1;
          int ext_start_bin = roi->get_ext_start_bin();
          int ext_end_bin = roi->get_ext_end_bin();
          if (ext_start_bin > roi1->get_start_bin()) ext_start_bin = roi1->get_start_bin();
          if (ext_end_bin < roi1->get_end_bin()) ext_end_bin = roi1->get_end_bin();
          roi->set_ext_start_bin(ext_start_bin);
          roi->set_ext_end_bin(ext_end_bin);
          //    std::cout << roi->get_chid() << " " << roi1->get_chid() << " " << roi->get_start_bin() << " " << roi->get_end_bin() << " " << roi1->get_start_bin() << " " << roi1->get_end_bin() << std::endl;
        }
        //std::cout << chid << " " << roi->get_start_bin() << " " << roi->get_end_bin() << std::endl;
      }

    }
    if (flag){
      SignalROI *prev_roi = 0;
      for (auto it = rois_w_tight.at(chid).begin(); it!= rois_w_tight.at(chid).end();it++){
        SignalROI *roi =  *it;
        if (prev_roi!=0){
          if (prev_roi->get_ext_end_bin() > roi->get_ext_start_bin()){
            prev_roi->set_ext_end_bin(int(prev_roi->get_end_bin() * 0.5 + roi->get_start_bin()*0.5));
            roi->set_ext_start_bin(int(prev_roi->get_end_bin() * 0.5 + roi->get_start_bin()*0.5));
          }
        }
        prev_roi = roi;
      }
    } 
      // for (auto it = rois_w_tight.at(chid).begin(); it!= rois_w_tight.at(chid).end();it++){
      //   SignalROI *roi =  *it;
      //   std::cout << chid << " " << roi->get_ext_start_bin() << " " << roi->get_ext_end_bin() << std::endl;
      // }
  }
}




// Local Variables:
// mode: c++
// c-basic-offset: 2
// End: