histgen.cxx

Go to the documentation of this file.

//
// This program is to generate histograms for sampling of dE/dx simulation
// First the optimized binning of beta-gamma of hadron is read from an external file
// No optimization for muon and electron, just use equi-width for them
// Then save the dE/dx hit values to tree with equi-width of angle
// Finally, binning them with equi-frequency 
//
 
#include "TCanvas.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1F.h"
#include <sstream>
#include <iostream>
#include <string>
#include <cstring>
#include <vector>
#include <algorithm>
#include <functional>
#include <cmath>
using namespace std;
// binning information for beta-gamma is read from file "binning.h", only modify it if necessary 
#include "binning.h"
 
const int na = 10;  // #angle
const int nc = 2;   // #charge
 
const double muon_step  = (double)(bg_muon_max - bg_muon_min) / (double)bg_muon_bin;
const double electron_step = (double)(bg_electron_max - bg_electron_min) / (double)bg_electron_bin;
 
static int histNum = 0;
 
int find_bgidx(int type, double bgg, double bg_min, double bg_step, int nbin);
double find_boundary(int type, int idx, double bg_min, double bg_step);
void t_to_h(TTree *t, TH1F *h);
void histgen(string, int, TObjArray&, vector<double>&);
void histps(string);
void hist_update(); // update histograms to the form suitable for database
const bool ps_flag(false);
const bool m_debug(false);
TFile *f_temp = new TFile("temp.root", "recreate");
 
int main(){
  std::cout << "in main()" << std::endl;
  //call fuction histgen to generate histograms
  TObjArray histCol;
  vector<double> bg;
  // now only three types, 0 for hadron, 3 for muon, 4 for electron
  std::cout << "before histgen" << std::endl;
  histgen("pre_data.root", 0, histCol, bg);
  std::cout << "hadron hists has been generated" << std::endl;
 
  histgen("pre_data.root", 3, histCol, bg);
  std::cout << "muon hists has been generated" << std::endl;
 
  histgen("pre_data.root", 4, histCol, bg);
  std::cout << "electron hists has been generated" << std::endl;
 
  //define file and write to it
  string filename = "DedxSampleHists.root";
  TFile* f = new TFile(filename.c_str(), "recreate");
 
  int total = histCol.GetLast() + 1;
  int bg_num = bg.size();
  if(m_debug) cout << "total: " << total << "  bg num = " << bg_num << endl;
  double betagamma[3000];
  for (unsigned int i = 0; i < bg.size(); i++) betagamma[i] = bg[i];
 
  TTree* t = new TTree("bin", "bin information");
  t->Branch("totalNum", &total, "total/I");
  t->Branch("betagammaBounds", &bg_num, "bg_num/I");
  t->Branch("betagamma", betagamma, "betagamma[bg_num]/D");
  t->Fill();
  t->Write();
    
  std::cout << "begin write hists" << std::endl;
  for (int i=0; i<total; i++){
    if(i%100==0) cout << "i: " << i << " histogram is written now" << endl;
    ((TH1F*)histCol[i])->Write();
  }
  std::cout << "write hists finish" << std::endl;
 
  f->Close();
  f_temp->Close();
  if(ps_flag) histps(filename);
  hist_update();
  return 0;
}
 
 
// find beta-gamma index by beta-gamma value
int find_bgidx(int type, double bgg, double bg_min, double bg_step, int nbin){
  if(type==0){
    for(int i=0; i<=nbin; i++){
      if(bgg<vec_bg_hadron[i+1]) return i;
    }
    cout << "not found correct bgidx, something is wrong!" << endl;
    return -1;
  }
  else return((int)(bgg-bg_min)/bg_step);
}
 
// find beta-gamma boundary by index 
double find_boundary(int type, int idx, double bg_min, double bg_step){
  if(type==0) return vec_bg_hadron[idx];
  else return(bg_min + bg_step*idx);
}
 
// generate histograms from tree with equi-frequency method
void t_to_h(TTree *t, TH1F *h){
  vector<double> vec_dedx(0);
  double m_dedx;
  t->SetBranchAddress("dedx_hit", &m_dedx);
  for(int i=0; i<t->GetEntries(); i++){
    t->GetEntry(i);
    vec_dedx.push_back(m_dedx);
  }
  sort(vec_dedx.begin(), vec_dedx.end());
  double points[101];
  int step = t->GetEntries()/100;
  int rad =  t->GetEntries()%100;
  int count = 0;
  int sum = 1;
  points[0] = vec_dedx.front();
  if(m_debug) cout << "points[0] " << points[0] << endl;
  for(unsigned int i=0; i<vec_dedx.size(); i++){
    count ++;
    if(count>=step){
      points[sum] = (vec_dedx[i]+vec_dedx[i+1])/2;
      sum ++;
      count = 0;
    }
    if(m_debug) cout << "i= " << i <<  "  rad " << rad << "  step " << step << "  count " << count << "  sum " << sum-1 << "  point " << points[sum-1] << endl;
    if(sum==100) break;
  }
  points[100] = vec_dedx.back();
  // if(m_debug) cout << "points[100] " << points[100] << endl;
  double points_mod[111];
  for(int j=0; j<111; j++){
    if(j<11) points_mod[j] = points[1]/11.*j;
    else points_mod[j] = points[j-10];
    if(m_debug) cout << "j= " << j <<  " modified  point " << points_mod[j] << endl;
  }
  h->SetBins(110, points_mod);
  t->Project(h->GetName(), "dedx_hit");
  t->Delete();
}
 
void histgen(string filename, int type, TObjArray& hist, vector<double>& bg){
  float exraw[100];
  float bgg, dEdx_meas, costheta, charge;
  int nbin, ndedxhit;
  double bg_min, bg_max, bg_step;
 
  TFile* f = new TFile(filename.c_str());
  TTree* trk = (TTree*)f->Get("n103");
  trk->SetBranchAddress("dEdx_meas",  &dEdx_meas);
  trk->SetBranchAddress("bg",     &bgg);
  trk->SetBranchAddress("costheta", &costheta);
  trk->SetBranchAddress("charge",    &charge);
  trk->SetBranchAddress("dEdx_hit",     exraw);
  trk->SetBranchAddress("nhits", &ndedxhit);
 
  //initialize variables
  switch (type) {
  case 0:
    nbin  = bg_hadron_bin;
    bg_min = bg_hadron_min;
    bg_max = bg_hadron_max;
    break;
  case 3:
    nbin  = bg_muon_bin;
    bg_min = bg_muon_min;
    bg_max = bg_muon_max;
    break;
  case 4:
    nbin  = bg_electron_bin;
    bg_min = bg_electron_min;
    bg_max = bg_electron_max;
    break;
  default:
    std::cout << "Error: no such sample type!" << std::endl;
  }
 
  //define trees and histograms
  f_temp->cd();
  TH1F** h = new TH1F*[nbin * 2 * na];
  TTree** t = new TTree*[nbin * 2 * na];
  stringstream ss, sst;
  double m_dedx_hit;
  for (int i = 0; i < nbin * 2 * na; i++){ // nbin is different for different type of particle
    ss.str(""); sst.str("");
    ss << "h" << histNum + i; sst << "t" << histNum + i;
    h[i] = new TH1F(ss.str().c_str(), ss.str().c_str(), 1, 0, 1);
    f_temp->cd();
    t[i] = new TTree(sst.str().c_str(), sst.str().c_str());
    t[i]->Branch("dedx_hit", &m_dedx_hit, "dedx_hit/D");
  }
  histNum += nbin*2*na; // the number of total histograms is a sum of No. of hists of each particle type
 
  //loop every hit of n103
  int bgidx(-1), angleidx(-1), chargeidx(-1), treeidx(-1);
    
  for (int i = 0; i < trk->GetEntries(); i++){
    trk->GetEntry(i);
    if(bgg<bg_min || bgg>bg_max) continue;
    //    if(m_debug) cout << "i " << i << "  dedx " << dEdx_meas << "  bgg " << bgg  << "  costheta " << costheta << "  charge " << charge << "  nhits " << ndedxhit << endl;
 
    // betagamma bin index
    bg_step = (bg_max-bg_min)/nbin;
    bgidx = find_bgidx(type, bgg, bg_min, bg_step, nbin);
    if(bgidx<0 || bgidx>=nbin) continue;
 
    // costheta bin index
    angleidx = (int)(fabs(costheta) * na / 0.93);
    if(angleidx<0 || angleidx>=na) continue;
 
    // charge bin index
    if (charge == 0) continue;
    else chargeidx = (charge > 0) ? 0 : 1;
    // if(m_debug) cout << "bgidx " << bgidx << "  angleidx " << angleidx << "  chargeidx " << chargeidx <<endl;
    
    // fill tree
    treeidx = -1;
    for (int j = 0; j <ndedxhit; j++) {
      //      for(int idx_charge = chargeidx; idx_charge< chargeidx+2; idx_charge++){
      treeidx = bgidx*na*2 + angleidx*2 + chargeidx;
      m_dedx_hit = exraw[j];
      //      if(m_debug) cout << "j " << j << "dedx_hit " << exraw[j] << endl;
      if(m_dedx_hit>0 && m_dedx_hit<500000) t[treeidx]->Fill(); // a very loose cut on dE/dx value
      //} 
    }// end of loop in all hits in one track
  }// end of loop in all tracks in one type of particle
 
  //push the histograms and betagammas to the vector
  for (int i = 0; i < nbin * na * nc; i++){
    if(t[i]->GetEntries()<1000) cout << "i " << i << "  no." << t[i]->GetEntries() << endl;
    t_to_h(t[i], h[i]);
    hist.Add(h[i]);
  }
  for (int i = 0; i < nbin; i++){
    bg.push_back(find_boundary(type, i, bg_min, bg_step));
    bg.push_back(find_boundary(type, i+1, bg_min, bg_step));
  }
}// end of histgen() function
 
 
void histps(string name){
  TFile f(name.c_str());
  TTree* t = (TTree*)f.Get("bin");
  Int_t total, bgb;
  t->SetBranchAddress("totalNum", &total);
  t->SetBranchAddress("betagammaBounds", &bgb);
  t->GetEntry(0);
  TH1F* h;
  stringstream ss;
  TCanvas c1("c1", "canvas", 800, 600);
  c1.Print((name + ".ps[").c_str());
  for (Int_t i = 0; i < total; i++) {
    if(i%100==0) cout << "i: " << i << " ps is written now" << endl;
    ss.str("");
    ss << "h" << i;
    h = (TH1F*)f.Get(ss.str().c_str());
    h->Draw();
    c1.Print((name + ".ps").c_str());
  }
  c1.Print((name + ".ps]").c_str());
}
 
 
void hist_update(){
  char hname[200];
  TH1F *h1 = new TH1F();
  TTree *tree = new TTree("TH1F_Col", "TH1F_Col");
  tree->Branch("TH1F_Col", "TH1F", &h1, 3200, 0);
  TFile* f1 = new TFile("DedxSampleHists.root", "UPDATE");
  TTree* bin = (TTree*)f1->Get("bin");
  Int_t totNum;
  bin->SetBranchAddress("totalNum", &totNum);
  bin->GetEntry(0);
 
  for (Int_t i = 0; i < totNum; i++){
    if(i >= 10000) sprintf(hname, "h%05d", i);
    if(i >= 1000)  sprintf(hname, "h%04d", i);
    if(i >= 100)   sprintf(hname, "h%03d", i);
    if(i >= 10)    sprintf(hname, "h%02d", i);
    if(i < 10)     sprintf(hname, "h%01d", i);
    h1 = (TH1F*)f1->Get(hname);
    tree->Fill();
  }
  tree->Write();
  f1->Close();
}