#include Double_t MeanInverse_func (Double_t *x, Double_t *par){ // returns average value of <1 + scale*(E1+E2)/E1*E2>, function of E1, average over E2. // constants Double_t Ethresh = par[0]; Double_t paLab0 = par[1]; Double_t scale=par[2]; Double_t E1 = x[0]; Double_t f = 0; if (E1 <= 0) return f; if (E1 > paLab0/2) { f = log( (paLab0-E1)/Ethresh) + (paLab0 - E1 - Ethresh)/E1; f = f / (paLab0-E1-Ethresh); } else { f = log(E1/Ethresh) + (E1-Ethresh)/E1; f = f / (E1 - Ethresh); } // cout << "E1=" << E1 << " f=" << f << endl; // f *= scale; f = scale*f<1? 1 - 0.5*scale*f: 0; // use 1/2 as approximation for square root return f; } Double_t Ecorr_func (Double_t *x, Double_t *par){ // parameterization of energy correction function: E_corr = E/Ecoor; // constants Double_t q0 = par[0]; Double_t q1 = par[1]; Double_t q2 = par[2]; Double_t q3 = par[3]; Double_t q4 = par[4]; Double_t q5 = par[5]; Double_t q6 = par[6]; Double_t q7 = par[7]; Double_t q8 = par[8]; Double_t Eg = x[0]; Double_t f = 0; f = Eg>0? q0 - q1*exp(q2*Eg+q3) - q4/(q5 + q6*exp(-q7*Eg+q8)) : 0; // cout << " Ecorr parameters=" << q0 << " " << q1 << " " << q2 << " " << q3 << " " << q4 << " " << q5 << " " << q6 << " " << q7 << " " << q8 << endl; // cout << "Ecorr_func: Eg=" << Eg << " f=" << f << endl; return f; } Double_t leakage_func (Double_t *x, Double_t *par){ // leakage parameterization based on BCAL NIM paper Fig. 10 // leakage as a fraction // constants Double_t k1 = par[0]; Double_t k2 = par[1]; Double_t k3 = par[2]; // Double_t ktheta = par[3]; Double_t ktheta = 90.; // eliminate theta dependence. Double_t k4 = par[4]; Double_t k5 = par[5]; Double_t Eg = x[0]; Double_t f1=0; Double_t f2=0; Double_t factor=4; // decrease angle dependence; nominal factor=1 Double_t lambda=factor*25; Double_t k1m = k1*exp(-(ktheta-12)/lambda); lambda=50*factor; Double_t factor2=4; Double_t k2m = factor2*k2*exp(-(ktheta-12)/lambda); lambda = 30*factor; Double_t k3m = k3*exp(-(ktheta-12)/lambda); // Front leakage f1 = Eg>0? k1m + k2m*exp(-Eg/k3m): 0; // return zero for negative Eg f1 = Eg>0? 0.2/Eg: 0; // fixed leakage lambda = 30; Double_t k4m = k4*exp(-(90-ktheta)/lambda); Double_t k5m = k5; // Back leakage // f2 = Eg>0? k4m*(1-exp(-Eg/k5m)) : 0; // constant at high energy f2 = Eg>0? k4m*sqrt(Eg/k5m) : 0; f2 = 0; // no leakage return 1- f1+f2; // return E1_measured/E1 = E1 *(1-L)/E1 = 1-L; } void plot_E1dep(void) { // plot the E1 dependence of the pi0 mass ratio for various scenarios // gStyle->SetPalette(1,0); gStyle->SetOptStat(111111); gStyle->SetOptFit(111111); gStyle->SetPadRightMargin(0.025); gStyle->SetPadLeftMargin(0.2); gStyle->SetPadBottomMargin(0.15); TString filename = "plot_E1dep"; // define function Double_t paLab0=8; Double_t Ethresh=0.5; Double_t k1=0.06; Double_t k2=0.1; // Double_t k3=0.5; // nominal value Double_t k3=4; // increase energy dependence Double_t ktheta=90; /* Double_t k4=0.03; // constants for exponential fit Double_t k5=0.3;*/ Double_t k4=0.03; // constants for exponential fit Double_t k5=1; Int_t npar = 6; Double_t xmin=1; Double_t xmax=7; TF1 *leakage = new TF1("leakage",leakage_func,xmin,xmax,npar); leakage->SetParameters(k1,k2,k3,ktheta,k4,k5); leakage->SetParNames("k1","k2","k3","ktheta","k4","k5"); Double_t fudge=1; TF1 *inverse = new TF1("inverse","[2]*(1 - 0.5*[0]*(x+[1])/(x*[1]))",xmin,xmax); inverse->SetParameters(0.2,1,fudge); inverse->SetParNames("leak","E2","fudge"); // fits to leakage fraction ~ 0.2/E. These should match at f=4. TF1 *e1 = new TF1("e1","[2]* (log( ([1]-x)/[0]) + ([1] - x - [0])/x)",1,4); TF1 *e2 = new TF1("e2","[2]* (1 + log(x/[0]) - 2*[0]/[1])",4,7); Double_t scale=1.2/9.; e1->SetParameters(Ethresh,paLab0,scale); e1->SetParNames("Ethresh","paLab0","scale"); e2->SetParameters(Ethresh,paLab0,scale); e2->SetParNames("Ethresh","paLab0","scale"); npar=3; TF1 *MeanInverse = new TF1("MeanInverse",MeanInverse_func,xmin,xmax,npar); scale=0.2; MeanInverse->SetParameters(Ethresh,paLab0,scale); MeanInverse->SetParNames("Ethresh","paLab0","scale"); TCanvas *c0 = new TCanvas("c0", "c0",200,10,1000,1000); c0->Divide(1,1); Double_t ymin=0.5; Double_t ymax=1.1; TString text; /*c0->cd(1); ktheta = 12.; text.Form("Leak=0.2/E, E_{#pi}=%.1f GeV, Ethresh=%.1f GeV",paLab0,Ethresh); TH2D *boundaries1 = new TH2D("boundaries1", "",xmax-xmin, xmin, xmax, ymax-ymin, ymin, ymax); boundaries1->SetTitle(text); boundaries1->SetTitleSize(0.07,"T"); boundaries1->SetStats(0); boundaries1->SetYTitle("m/m_{#pi}"); boundaries1->SetXTitle("E1 (GeV)"); boundaries1->GetYaxis()->SetLabelSize(0.07); boundaries1->GetYaxis()->SetTitleSize(0.07); boundaries1->GetYaxis()->SetTitleOffset(1.3); boundaries1->GetXaxis()->SetLabelSize(0.07); boundaries1->GetXaxis()->SetTitleSize(0.07); boundaries1->GetXaxis()->SetTitleOffset(1.0); boundaries1->GetXaxis()->SetNdivisions(510); boundaries1->Draw();*/ TF1 *leakage_draw = new TF1("leakage_draw",leakage_func,xmin,xmax,npar); leakage_draw->SetParameters(k1,k2,k3,ktheta,k4,k5); leakage_draw->SetParNames("k1","k2","k3","ktheta","k4","k5"); // leakage_draw->Draw("same"); /*c0->cd(2); e1->Draw(); e2->Draw("same");*/ /*c0->cd(3); gPad->SetGridx(); gPad->SetGridy(); boundaries->Draw(); e1->Draw("same"); e2->Draw("same"); inverse->Draw("same");*/ ymin=0.5; ymax=1.1; // TString text; text.Form("Leak=0.2/E, E_{#pi}=%.1f GeV, Ethresh=%.1f GeV",paLab0,Ethresh); TH2D *boundaries = new TH2D("boundaries", "",xmax-xmin, xmin, xmax, ymax-ymin, ymin, ymax); boundaries->SetTitle(text); boundaries->SetTitleSize(0.07,"T"); boundaries->SetStats(0); boundaries->SetYTitle("m/m_{#pi}"); boundaries->SetXTitle("E1 (GeV)"); boundaries->GetYaxis()->SetLabelSize(0.07); boundaries->GetYaxis()->SetTitleSize(0.07); boundaries->GetYaxis()->SetTitleOffset(1.3); boundaries->GetXaxis()->SetLabelSize(0.07); boundaries->GetXaxis()->SetTitleSize(0.07); boundaries->GetXaxis()->SetTitleOffset(1.0); boundaries->GetXaxis()->SetNdivisions(510); c0->cd(1); gPad->SetGridx(); gPad->SetGridy(); boundaries->Draw(); MeanInverse->SetLineColor(1); MeanInverse->Draw("same"); TF1 *inverse1 = inverse->DrawCopy("same"); fudge = 1; inverse1->SetParameters(0.2,1,fudge); inverse1->SetRange(1,paLab0-1); inverse1->SetLineStyle(1); inverse1->SetLineColor(4); TF1 *inverse2 = inverse->DrawCopy("same"); fudge = 1; inverse2->SetParameters(0.2,2,fudge); inverse2->SetRange(2,paLab0-2); inverse2->SetLineStyle(2); inverse2->SetLineColor(4); TF1 *inverse3 = inverse->DrawCopy("same"); fudge = 1; inverse3->SetParameters(0.2,3,fudge); inverse3->SetRange(3,paLab0-3); inverse3->SetLineStyle(3); inverse3->SetLineColor(4); leakage_draw->Draw("same"); leakage_draw->SetLineColor(2); TLegend *leg = new TLegend(0.3,0.2,0.6,0.4); leg->AddEntry(leakage,"1 - Leakage","l"); leg->AddEntry(inverse1,"Ratio (E2=1 GeV)","l"); leg->AddEntry(inverse2,"Ratio (E2=2 GeV)","l"); leg->AddEntry(inverse3,"Ratio (E2=3 GeV)","l"); leg->AddEntry(MeanInverse,"Ratio ","l"); leg->Draw(); /*c2->cd(6); h1_DeltaCos->SetTitle(text); // h1_E2measVSE1meas->GetXaxis()->SetRangeUser(xmin,xmax); // h1_DeltaCos->GetYaxis()->SetRangeUser(ymin,ymax); h1_DeltaCos->GetXaxis()->SetTitleSize(0.05); h1_DeltaCos->GetYaxis()->SetTitleSize(0.05); h1_DeltaCos->GetXaxis()->SetTitle("CosTheCM_meas - CosTheCM"); h1_DeltaCos->GetYaxis()->SetTitle("Events"); h1_DeltaCos->SetLineColor(1); h1_DeltaCos->Fit("gaus"); h1_DeltaCos->Draw();*/ // TString outfile; // outfile.Form(filename+"_%.0f_%.0f_%.0f_%.0f_%d",1000*paLab0,1000*Ethresh,1000*BCALa,1000*BCALb,iter); c0->SaveAs(filename+".pdf"); }