prepare_data.cxx

Go to the documentation of this file.

/* 1) This program is to collect all types of particles in an uniform style (dE/dx, angle, beta-gamma)
    Output the data to a preliminary.root file
 2) Also do a coarse region devision, the results should be output to a preliminary.txt file 
    with other additional information
*/
#include "TCanvas.h"
#include "TFile.h"
//#include "TKey.h"
#include "TTree.h"
#include "TH1F.h"
//#include <TROOT.h>
#include <sstream>
#include <iostream>
#include <string>
#include <cstring>
#include <vector>
#include <cmath>
#include "TCut.h"
 
using namespace std;
 
#define mass_e          0.511e-3
#define mass_mu         105.658e-3  
#define mass_pi         139.570e-3
#define mass_k          493.677e-3
#define mass_p          938.272e-3 
 
float  ptrk, dEdx_meas, costheta, charge, bg, mass; // bg means beta-gamma
int nhits;
float dEdx_hit[100];
void collect(TString filename, int type, TTree *tree);
void setregion();
 
int main(){
  TFile *f = new TFile("pre_data.root","recreate");
  TTree *t = new TTree("n103","predata");
  t->Branch("dEdx_meas", &dEdx_meas, "dEdx_meas/f");
  t->Branch("costheta", &costheta, "costheta/f");
  t->Branch("charge", &charge, "charge/f");
  t->Branch("bg", &bg, "bg/f");
  t->Branch("nhits", &nhits, "nhits/I");
  t->Branch("dEdx_hit", dEdx_hit, "dEdx_hit[nhits]/f");
 
  // collect all data here
  cout << "begin to collect data" << endl;
  collect("../hadron_track/proton/proton.root", 0, t);
  cout << "proton collected" << endl;
  collect("../hadron_track/kaon/kaon.root", 1, t);
  cout << "kaon collected" << endl;
  collect("../hadron_track/pion/pion.root", 2, t);
  cout << "pion collected" << endl;
  collect("../hadron_track/muon/muon.root", 3, t);
  cout << "muon collected" << endl;
  collect("../hadron_track/electron/electron.root", 4, t);
  cout << "electron collected" << endl;
 
  f->cd();
  t->Write();
  f->Close();
 
  //  setregion();
 
  return 0;
}
 
void collect(TString filename, int type, TTree *t){
  TFile *ff = new TFile(filename);
  TTree* trk = (TTree*)ff->Get("n103");
  trk->SetBranchAddress("dEdx_meas",  &dEdx_meas);
  trk->SetBranchAddress("ptrk",     &ptrk);
  trk->SetBranchAddress("costheta", &costheta);
  trk->SetBranchAddress("charge",    &charge);
  trk->SetBranchAddress("ndedxhit",  &nhits);
  trk->SetBranchAddress("dEdx_hit",   dEdx_hit);
 
  //initialize variables
  switch (type) {
  case 0:
    mass  = mass_p;
    break;
  case 1:
    mass  = mass_k;
    break;
  case 2:
    mass  = mass_pi;
    break;
  case 3:
    mass  = mass_mu;
    break;
  case 4:
    mass  = mass_e;
    break;
  default:
    cout << "Error: no such sample type!" << endl;
  }
 
  int m_ntot = trk->GetEntries();    
  // if(m_ntot>500000) m_ntot = 500000;
  for (int i = 0; i < m_ntot; i++){
    trk->GetEntry(i);
    if (type == 0 && dEdx_meas < 500) continue; // type 0 is proton!
    if (type == 4 && dEdx_meas > 1200) continue; // type 4 is electron!
    if(ptrk <= 0 || ptrk > 2.3) continue;
    if(fabs(costheta) >= 0.93) continue;
    bg = ptrk/mass;
    t->Fill();
  }
  ff->Close();
}// end of collect()
 
void setregion(){ // use greedy algorithm, step by step to check
  TFile *f = new TFile("pre_data.root");
  TTree *t = (TTree*)f->Get("n103");
  bool small_flag(true);
  double step(0.1); // for beta-gamma
  double current(0.0); // current position in beta-gamma
  int accu(0); // accumulated events
  double stop(4600);
  int threshold(1000);
  while(current<stop){
    cout << "now: " << current << endl;
    if(small_flag){
      accu = t->GetEntries(Form("bg>=%f && bg<%f", current, current+step));
      current += step;
      if(accu > threshold){
    cout << "point: " << current << "  events: " << accu << endl;
    small_flag = false;
      }
    }// end of if(small_flag)
    else{
      accu = t->GetEntries(Form("bg>=%f && bg<%f", current, current+step));
      current += step;
      if(accu < threshold){
    cout << "point: " << current << "  events: " << accu << endl;
    small_flag = true;
      }
    }// end of else(small_flag)
  }// end of while
  f->Close();
}