00001 #include <iostream>
00002 #include <fstream>
00003 using namespace std;
00004
00005 #include "TFile.h"
00006 #include "TMath.h"
00007 #include "TH1.h"
00008 #include "TF1.h"
00009 #include "TPostScript.h"
00010 #include "TCanvas.h"
00011 #include "TStyle.h"
00012 #include "TLatex.h"
00013 #include "TString.h"
00014 #include "TLine.h"
00015
00016 void drawTower(TH1* h, TF1* f, int id);
00017
00018 void mip_ring_draw(const char* infile="mipskimfile.root", const char* postscript="mipring.ps", const char* outname="2008mipring.root")
00019 {
00020 const int nrings = 40;
00021 TFile input(infile);
00022 char name[100];
00023 TH1* ring_histo[nrings];
00024 TH1F* ringMIP = new TH1F("ringMIP","",nrings,-1.0,1.0);
00025 ringMIP->SetXTitle("#eta");
00026 ringMIP->SetYTitle("MIP peak (ADC)");
00027 for(int i = 0; i < nrings; i++)
00028 {
00029 sprintf(name,"ring_histo_%i",i+1);
00030 ring_histo[i] = (TH1F*)input.Get(name);
00031 }
00032
00033
00034 TPostScript* ps = new TPostScript(postscript);
00035 TCanvas* c = new TCanvas("c","",100,100,600.,800.);
00036
00037
00038
00039
00040 int pad = 16;
00041 TF1 *gaussian_fit[nrings], *landau_fit[nrings];
00042 for(int i=0; i<nrings; i++){
00043 if(pad%15 == 1){
00044 c->Update();
00045 ps->NewPage();
00046 c->Clear();
00047 c->Divide(3,5);
00048 pad = 1;
00049 }
00050 c->cd(pad);
00051
00052 cout<<"fitting ring "<<i+1<<" of "<<nrings<<endl;
00053
00054 sprintf(name,"fit_%i",i+1);
00055
00056 gaussian_fit[i] = new TF1(name,"gaus",7.,45.);
00057 landau_fit[i] = new TF1(name,"landau",7.,200.);
00058
00059 gaussian_fit[i]->SetParameter(1,20.);
00060 gaussian_fit[i]->SetParameter(2,5.);
00061
00062 landau_fit[i]->SetParameter(1,17.);
00063 landau_fit[i]->SetParameter(2,3.);
00064
00065 gaussian_fit[i]->SetLineColor(kGreen);
00066 gaussian_fit[i]->SetLineWidth(0.6);
00067
00068 landau_fit[i]->SetLineColor(kYellow);
00069 landau_fit[i]->SetLineWidth(0.6);
00070
00071 ring_histo[i]->Fit(gaussian_fit[i],"rq");
00072 ring_histo[i]->GetXaxis()->SetRangeUser(3.0,50.0);
00073
00074 drawTower(ring_histo[i],gaussian_fit[i],i);
00075
00076 double histogram_top = ring_histo[i]->GetBinContent(ring_histo[i]->GetMaximumBin());
00077
00078 double gaussian_mean = gaussian_fit[i]->GetParameter(1);
00079 TLine *gaussian_peak = new TLine(gaussian_mean,0.,gaussian_mean,histogram_top+15);
00080 gaussian_peak->SetLineColor(kGreen);
00081 gaussian_peak->SetLineWidth(2.0);
00082 gaussian_peak->Draw("same");
00083
00084 pad++;
00085 float eta = (float)((i - 20) * 2 + 1)/40;
00086 float theta = 2*TMath::ATan(TMath::Exp(-eta));
00087 float mipenergy = 0.264 * (1+0.056*eta*eta)/TMath::Sin(theta);
00088 cout<<mipenergy/gaussian_mean<<", ";
00089 ringMIP->Fill(eta,gaussian_mean);
00090 ringMIP->SetBinError(i+1,gaussian_fit[i]->GetParError(1));
00091
00092 }
00093 cout<<endl;
00094 c->Update();
00095 ps->NewPage();
00096 c->Clear();
00097 c->Divide(1,2);
00098 c->cd(1);
00099 ringMIP->SetMarkerStyle(20);
00100 ringMIP->GetYaxis()->SetRangeUser(0,35);
00101 gStyle->SetOptStat(0);
00102 ringMIP->Draw("ep");
00103
00104 ps->Close();
00105
00106 TFile outfile(outname,"RECREATE");
00107 for(int i = 0; i < nrings; i++)
00108 {
00109 ring_histo[i]->Write();
00110 }
00111 ringMIP->Write();
00112 outfile.Close();
00113 }
00114
00115 void drawTower(TH1* h, TF1* f, int id){
00116
00117
00118
00119
00120
00121
00122 int xLatexBin = 20;
00123
00124
00125 h->SetXTitle("ADC");
00126 h->Draw();
00127
00128
00129 char line_name[50];
00130 sprintf(line_name,"mip_peak_%i",id);
00131
00132
00133
00134
00135
00136
00137 char tower_title[100];
00138 float eta = (float)((id - 20) * 2 + 1)/40;
00139 sprintf(tower_title,"eta = %f",eta);
00140 TLatex title_latex;
00141 title_latex.SetTextSize(0.15);
00142
00143
00144 title_latex.DrawTextNDC(0.13,0.78,tower_title);
00145
00146
00147 }