binning.cxx

Go to the documentation of this file.

//
// This program is to optimize the binning, a mixture method is used
// angle: equi-width, dedx: equi-frequency, beta-gamma: maximize information entropy
// for information entropy method, firstly we should classified data into different groups with respect to beta-gamma
//
 
#include "TCanvas.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1F.h"
#include "TMath.h"
#include "TObjectTable.h"
#include <sstream>
#include <iostream>
#include <fstream>
#include <string>
#include <cstring>
#include <vector>
#include <cmath>
 
using namespace std;
 
TTree* t = new TTree("t", "tree to save original data");
double m_left(0), m_right(0);
int m_bins(0);
typedef std::vector<double> vdouble;
bool m_debug(true);
vdouble vec_bg(0);
vdouble v(0); // a temp vector
void f_opt(double left, double right, int bins); // optimize binning, find beta-gamma points
double f_cals(double left, double right, double step); // calculate entropy  
void f_gen_hists(){};
void print_vec(vdouble v);
void print_vec_verso(vdouble v);
const double test_times(2);
double total_n(0); // total nhits in one region
TTree *t_save = new TTree("t_save", "tree to save entropy"); // used to save calculated entropy
double save_left(0), save_right(0), save_pos(0), save_s(0);
int tot_num_group; // number of groups
double x_group[2000]; // save the boundary of groups
double hits_group[2000]; // save hits in a group
 
int main(int argc, char **argv){
  //  gObjectTable->Print();
  if(argc<2){ cout << "please append config file, ex. binning.exe config.txt" << endl; return 0;}
  TFile *f = new TFile("pre_data.root");
  t = (TTree*)f->Get("n103");
  ifstream config_file(argv[1]);
  if(!config_file){
    cout << "cannot read config file correctly!" << endl;
    return 1;
  }
  TFile *f_group = new TFile("group_points.root");
  if(!f_group){
    cout << "cannot read group points file correctly!" << endl;
    return 1;
  }
  TTree *t_group = (TTree*)f_group->Get("t");
  double m_x_temp(0);
  t_group->SetBranchAddress("x", &m_x_temp);
  tot_num_group = t_group->GetEntries();
  for(int i=0; i<tot_num_group; i++) { t_group->GetEntry(i); x_group[i] = m_x_temp; }
  // t_group->Delete();
  f_group->Close();
  delete f_group;
  // prepare tree t_save
  t_save->Branch("save_left", &save_left, "save_left/D");
  t_save->Branch("save_right", &save_right, "save_right/D");
  t_save->Branch("save_pos", &save_pos, "save_pos/D");
  t_save->Branch("save_s", &save_s, "save_s/D");
  //
  while(m_bins!=-1){ // -1 is the end signal
    config_file >> m_left >> m_right >> m_bins >> m_debug;
    if(m_debug) cout << "config: " << m_left << " " << m_right << " " << m_bins << endl;
    f_opt(m_left, m_right, m_bins);
    cout <<"finial results for this time binning" << endl;
    print_vec_verso(vec_bg);
    f_gen_hists();
  }
  t->Delete();
  t_save->Delete();
  f->Close();
  delete f; 
  return 0;
}
 
 
void f_opt(double m_left, double m_right, int bins){
  if(m_debug) cout << "f_opt() is happily running!" << endl;
  int total(1); // number of bins
  vec_bg.clear();
  vec_bg.push_back(m_left);
  vec_bg.push_back(m_right);
  TTree *tree = new TTree("tree", "tree to save temp data");
  tree = t->CopyTree(Form("bg>%f && bg<%f", m_left, m_right));
  int m_hits;
  total_n = 0;
  tree->SetBranchAddress("nhits", &m_hits);
  int m_num_total(tree->GetEntries());
  for(int k=0; k<m_num_total; k++){
    tree->GetEntry(k);
    total_n += m_hits;
  }
  tree->Delete();
 
  while(total<bins){
    double temp_s(0), temp_pos(0);
    int temp_i(-1);
    for(unsigned int i=0; i<vec_bg.size(); i++){
      double pos_step((vec_bg[i+1]-vec_bg[i])/test_times);
      double ori_s(-1);
      for(double pos=vec_bg[i]+pos_step; pos<vec_bg[i+1]-0.00001; pos+=pos_step){
    if(ori_s<0) ori_s = f_cals(vec_bg[i], vec_bg[i+1], vec_bg[i+1]);// for the ori one, no mid point
    v.clear();
    v = vec_bg;
    v.insert(v.begin()+i+1, pos);
    double diff_s = ori_s - f_cals(vec_bg[i], vec_bg[i+1], pos); // information gain 
    if(diff_s > temp_s){
      temp_s = diff_s;
      temp_i = i;
      temp_pos = pos;
    }
    if(m_debug){
      cout << "i:" << i << "  pos:" << pos << "  diff_s:" << diff_s << "  temp_s:" << temp_s << endl;
      cout << "temp_v: " << endl;
      print_vec(v);
      cout << "vec_bg: " << endl;
      print_vec(vec_bg);
    }// end of print results
      }// end of loop in inert region
    }// end of loop in whole region
    vec_bg.insert(vec_bg.begin()+temp_i+1, temp_pos);
    cout << "============== vec_bg =============" <<endl;
    print_vec(vec_bg);
    cout << endl << endl;
    total ++;
  }// end of while()
}// end of f_opt()
 
 
void f_group(double left, double right, int &num, int &index){
  num = 0;
  index = -1;
  for(int i=1; i<tot_num_group; i++){
    if(num==0 && left<x_group[i]) { index = i-1; num = 1; }
    if(right<=x_group[i]){ num = i - index; break;}
  }
}
 
int f_group(double bg){
  for(int i=1; i<tot_num_group; i++){
    if(bg<x_group[i]) return i-1; 
  }
  return tot_num_group;
}
 
double f_cals(double left, double right, double pos){ 
  //  gObjectTable->Print();
  if(m_debug) cout << "f_cal() is happyily running!" << endl;
  for(int i=0; i<t_save->GetEntries(); i++){
    t_save->GetEntry(i);
    if(fabs(save_left-left)<0.0001 && fabs(save_right-right)<0.0001 && fabs(save_pos-pos)<0.0001) return save_s;
  } 
  double left_s(0), right_s(0), total_s(0);
  double left_ratio(0), right_ratio(0);
  TFile *file = new TFile("temp.root", "recreate");
  TTree *tree = t->CopyTree(Form("bg>%f && bg<%f", left, right));
  if(m_debug){ cout << "after open temp.root" << endl; gDirectory->pwd();  gDirectory->ls(); }
  double a_step(0.93*2/10);
  for(int i=0; i<2; i++){// only do test twice
    if(i==0) {m_left = left; m_right = pos;}
    else {m_left = pos; m_right = right;}
    int m_num_group(0);// how many groups will be covered by this bin
    int m_index_group(-1); // the index of left boundary of the groups
    f_group(m_left, m_right, m_num_group, m_index_group);
    if(m_debug) cout << "left " << m_left << " right " << m_right << " num of groups: " << m_num_group << "  boundary index: " << m_index_group << endl;
    for(int j=0; j<10; j++){// index of angle
      TTree *temp_tree = new TTree("tepm_tree", "temp tree");
      temp_tree = tree->CopyTree(Form("bg>%f && bg<%f && costheta>%f && costheta<%f", m_left, m_right, -0.93+j*a_step, -0.93+(j+1)*a_step));
      int nhits(0);
      Float_t m_bg(0);
      double tot_nhits(0);
      temp_tree->SetBranchAddress("nhits", &nhits);
      temp_tree->SetBranchAddress("bg", &m_bg);
      for(int m=0; m<tot_num_group; m++) hits_group[m]=0; 
      for(int k=0; k<temp_tree->GetEntries(); k++){// loop in track
    temp_tree->GetEntry(k);
    tot_nhits += nhits;
    hits_group[f_group(m_bg)] += nhits;
      }
      if(tot_nhits<2000){ file->Close(); return(100000.); } // if too small hits in a bin, bad (large) entropy
      else{
    for(int n=0; n<m_num_group; n++){ // loop in groups
      if(hits_group[m_index_group+n]==0) continue;
      // cout << "group_hits " << hits_group[m_index_group+n] << "  tot_nhits " << tot_nhits << endl;
      double p(hits_group[m_index_group+n]/tot_nhits); // the probability to find a hit of group_k in a bin
      if(p>1) cout << "m_index_group " << m_index_group << "  n " << n << " hits_group[] " <<  hits_group[m_index_group+n] << "  tot_nhits " << tot_nhits << endl;
      if(i==0) left_s += -p*TMath::Log2(p);
      if(i==1) right_s += -p*TMath::Log2(p);
    }// end of loop in groups
    if(i==0) left_ratio += tot_nhits;
    if(i==1) right_ratio += tot_nhits;
      }
      if(m_debug) cout << "m_left:" << m_left << "  m_right:" << m_right << "  j:" << j << "  left_s:" << left_s << "  right_s:" << right_s << endl;
      // temp_tree->Delete();
    }// end of loop angle
    if(fabs(right-pos)<0.00001) break;
  }// end of loop beta-gamma
  //  tree->Delete();
  file->Close();
  left_ratio = left_ratio/(left_ratio+right_ratio);
  right_ratio = 1-left_ratio;
  total_s = left_ratio*left_s + right_ratio*right_s;
  if(m_debug) cout << "left_ratio: " << left_ratio << "  right_ratio: " << right_ratio << "  total_s: " << total_s << endl;
  // save to the tree t_save
  save_left = left;
  save_right = right;
  save_pos = pos;
  save_s = total_s;
  t_save->Fill();
  return(total_s);
}
 
void print_vec(vdouble v){
  for(vdouble::iterator i=v.begin(); i<v.end(); i++) cout << *i << "  " ;
  cout << endl;
}
 
void print_vec_verso(vdouble v){
  cout << "vec[" << m_bins << "] = {" ;
  for(vdouble::iterator i=v.begin(); i<v.end(); i++) cout << *i << ", " ;
  cout << "};" << endl;
}