drawmom.cxx

Go to the documentation of this file.

#include "betagamma.cxx"
TString str_base = "../hadron_track/";
TString str_vec[5] = {"electron","muon","pion","kaon","proton"};
const double mom_max_vec[5] = {1.843, 2.0, 1.5, 1.8, 1.2};
double bg_max_vec[5];
const double masses[5] = {0.511e-3, 105.658e-3, 139.570e-3, 493.677e-3, 938.272e-3};
const double dedx_max_vec[5] = {1000, 1000, 1000, 1500, 7000};
double bg_matrix[5][2];
double mom_matrix[5][2];
const double y_matrix[5][2] = {{1.0,1.0},{2.0,2.0},{3.0,3.0},{4.0,4.0},{5.0,5.0}};
bool m_debug = false;
const TCut vec_cut[5]={"","","","","dEdx_meas>500"};
 
void drawmom(){
  for(int i=0; i<5; i++) bg_max_vec[i]=mom_max_vec[i]/masses[i];
  cout << "begin of drawmom()" << endl;
  TTree *t[5];
  TFile *f[5];
  TH1F *h[5];
  TGraph *g[5];
  TMultiGraph *mg = new TMultiGraph();
  TLatex *l[5];
  TCanvas *c1 = new TCanvas("momentum","momentum",10,10,1200,642);
  c1->Divide(3,2);
  // draw momentum histograms
  for(int i=0; i<5; i++){
    f[i] = new TFile(str_base + str_vec[i] + "/" + str_vec[i] + ".root");
    t[i] = (TTree*)f[i]->Get("n103");
    h[i] = new TH1F(str_vec[i], str_vec[i], (int)((mom_max_vec[i]-0.0)/0.001), 0.0, mom_max_vec[i]);
    c1->cd(i+1);
    t[i]->Project(str_vec[i], "ptrk", vec_cut[i], "", 100000);
    h[i]->Draw();
    mom_matrix[i][0] = h[i]->GetBinCenter(h[i]->FindFirstBinAbove(6.0));
    mom_matrix[i][1] = h[i]->GetBinCenter(h[i]->FindLastBinAbove(6.0));
    bg_matrix[i][0] = bg(i,mom_matrix[i][0]);
    bg_matrix[i][1] = bg(i,mom_matrix[i][1]);
    g[i] = new TGraph(2, bg_matrix[i], y_matrix[i]);
    mg->Add(g[i]);
    l[i] = new TLatex(bg_matrix[i][1], y_matrix[i][1], str_vec[i]);
  }
  c1->SaveAs("figs/mom.eps");
 
  // draw beta-gamma histograms
  delete c1;
  TCanvas *c1 = new TCanvas("#beta#gamma","#beta#gamma",10,10,1200,642);
  c1->Divide(3,2);
  for(int i=0; i<5; i++){
    delete h[i];
    h[i] = new TH1F(str_vec[i], str_vec[i], (int)((bg_max_vec[i]-0.0)/0.001), 0.0, bg_max_vec[i]);
    c1->cd(i+1);
    t[i]->Project(str_vec[i], Form("ptrk/%f",masses[i]), vec_cut[i], "", 100000);
    if(i) h[i]->Draw();
  }
  c1->SaveAs("figs/bg_hists.eps");
 
  TCanvas *c2 = new TCanvas("betagamma","betagamma",10,10,800,642);
  c2->SetLogx();
  mg->Draw("al");
  mg->GetXaxis()->SetTitle("#beta#gamma range");
  for(int j=0; j<5; j++) l[j]->Draw();
  c2->SaveAs("figs/bg.eps");
  // output information
  if(m_debug){
    for(int i=0; i<5; i++){
      cout << str_vec[i] << "  bins " << h[i]->GetNbinsX() << endl;
      cout << "mom  low: " << mom_matrix[i][0] << "  high: " << mom_matrix[i][1] << endl;
      cout << "bg  low: " << bg_matrix[i][0] << "  high: " << bg_matrix[i][1] << endl;
    }
  }
 
  TCanvas *c3 = new TCanvas("dedx vs momentum","dedx vs momentum",10,10,1200,642);
  c3->Divide(3,2);
  for(int i=0; i<5; i++){
    c3->cd(i+1);
    if(i==4) t[i]->Draw("dEdx_meas:ptrk", Form("ptrk<%f & dEdx_meas<%f & dEdx_meas>500", mom_max_vec[i], dedx_max_vec[i]), "", 100000);
    else t[i]->Draw("dEdx_meas:ptrk", Form("ptrk<%f & dEdx_meas<%f", mom_max_vec[i], dedx_max_vec[i]), "", 100000);
  }
  c3->SaveAs("figs/dedx_vs_mom.eps");
 
}
 
const double input_bg[5][2] = {{421.97,3187.82}, {10.98,115.6}, {0.74,10.98}, {0.64,0.74}, {0.26,0.64}};
void banch_rebg(){
  double output_mom[5][2];
  for(int i=0; i<5; i++){
    output_mom[i][0] = rebg(i,input_bg[i][0]);
    output_mom[i][1] = rebg(i,input_bg[i][1]);
    cout << str_vec[i] << endl;
   }
  for(int j=0; j<5; j++)  cout << "low " << output_mom[j][0] << "  high " << output_mom[j][1] << "  virtual bins for 1MeV width:" << (output_mom[j][1]-output_mom[j][0])/0.001 << endl << endl;
 
}