smd_offlmon/macro/plot_asym.C

0001 void plot_simple(TPad* pad, TH1* diff, TH1* sum, const char* which) {
0002     TF1* fit = new TF1("fitsimple", "[0] + [1]*sin(x - [2])", -3.141592, 3.141592);
0003     fit->SetParLimits(1, 0.001, 0.2);
0004     fit->SetParameter(1, 0.01);
0005     /* fit->SetParLimits(1, 3.041592, 3.241592); */
0006     fit->SetParLimits(2, 0.0, 2*3.141592);
0007     fit->SetParameter(2, 3.141592);
0008     fit->SetParName(0, "offset");
0009     fit->SetParName(1, "#epsilon");
0010     fit->SetParName(2, "#phi_{0}");
0011
0012     TH1F* h_asym = (TH1F*)diff->Clone("h_asym");
0013     h_asym->SetTitle(Form("%s SMD Simple Asymmetry", which));
0014     h_asym->SetYTitle("Raw Asymmetry");
0015     h_asym->Divide(diff, sum, 1.0, 1.0, "B"); // "B" option = binomial errors
0016
0017     pad->cd();
0018     gStyle->SetOptFit();
0019     gStyle->SetOptStat(0);
0020     h_asym->Draw();
0021     h_asym->Fit("fitsimple");
0022     float offset = fit->GetParameter(0);
0023     h_asym->GetYaxis()->SetRangeUser(offset-0.03, offset+0.03);
0024     TLine* l = new TLine(h_asym->GetXaxis()->GetXmin(), offset, h_asym->GetXaxis()->GetXmax(), offset);
0025     l->SetLineColor(kRed);
0026     l->SetLineStyle(kDashed);
0027     l->Draw();
0028     pad->Update();
0029 }
0030
0031 void plot_sqrt(TPad* pad, TH1* upleft, TH1* upright, TH1* downleft, TH1* downright, const char* which) {
0032     TF1* fit = new TF1("fitsqrt", "[0]*sin(x - [1])", -3.141592, 0.0);
0033     fit->SetParLimits(0, 0.001, 0.2);
0034     fit->SetParameter(0, 0.01);
0035     /* fit->SetParLimits(1, 3.041592, 3.241592); */
0036     fit->SetParLimits(1, 0.0, 2*3.141592);
0037     fit->SetParameter(1, 3.141592);
0038     fit->SetParName(0, "#epsilon");
0039     fit->SetParName(1, "#phi_{0}");
0040
0041     int nbins = upleft->GetNbinsX();
0042     TGraphErrors* sqrt_asym = new TGraphErrors(upleft);
0043     sqrt_asym->SetTitle(Form("%s SMD Square Root Asymmetry", which));
0044     sqrt_asym->GetXaxis()->SetTitle("#phi");
0045     sqrt_asym->GetYaxis()->SetTitle("Raw Asymmetry");
0046     for (int i=0; i<nbins; i++) {
0047     float n1, n2, n3, n4, e1, e2, e3, e4;
0048     n1 = upleft->GetBinContent(i+1);
0049     n2 = upright->GetBinContent(i+1);
0050     n3 = downleft->GetBinContent(i+1);
0051     n4 = downright->GetBinContent(i+1);
0052     e1 = upleft->GetBinError(i+1);
0053     e2 = upright->GetBinError(i+1);
0054     e3 = downleft->GetBinError(i+1);
0055     e4 = downright->GetBinError(i+1);
0056
0057     float asym = (sqrt(n1*n4)-sqrt(n3*n2))/(sqrt(n1*n4)+sqrt(n3*n2));
0058     float t1 = sqrt(n2*n3*n4/n1)*e1;
0059     float t2 = sqrt(n1*n3*n4/n2)*e2;
0060     float t3 = sqrt(n1*n2*n4/n3)*e3;
0061     float t4 = sqrt(n1*n2*n3/n4)*e4;
0062     float asym_err = (1/pow(sqrt(n1*n4)+sqrt(n3*n2),2))*sqrt(pow(t1,2)+pow(t2,2)+pow(t3,2)+pow(t4,2));
0063
0064     sqrt_asym->SetPointY(i, asym);
0065     sqrt_asym->SetPointError(i, sqrt_asym->GetErrorX(i), asym_err);
0066     }
0067     pad->cd();
0068     gStyle->SetOptFit();
0069     gStyle->SetOptStat(0);
0070     sqrt_asym->Draw("ap");
0071     sqrt_asym->Fit("fitsqrt");
0072     sqrt_asym->GetYaxis()->SetRangeUser(-0.03, 0.03);
0073     pad->Update();
0074 }
0075
0076 void plot_xy(TCanvas* c, TH2* northxy, TH2* southxy, TH1** northx, TH1** southx) {
0077     c->Clear();
0078     gStyle->SetOptStat(1);
0079     c->Divide(2, 2, 0.025, 0.025);
0080
0081     c->cd(1);
0082     /* gPad->SetLogz(); */
0083     northxy->Draw("colz");
0084     gPad->Update();
0085     TPaveStats* stats = (TPaveStats*)gPad->GetPrimitive("stats");
0086     stats->SetX1NDC(0.15);
0087     stats->SetX2NDC(0.295);
0088     stats->SetOptStat(1110);
0089
0090     c->cd(2);
0091     /* gPad->SetLogz(); */
0092     southxy->Draw("colz");
0093     gPad->Update();
0094     stats = (TPaveStats*)gPad->GetPrimitive("stats");
0095     stats->SetX1NDC(0.15);
0096     stats->SetX2NDC(0.295);
0097     stats->SetOptStat(1110);
0098
0099     c->cd(3);
0100     northx[0]->Sumw2();
0101     northx[0]->Scale(1.0/northx[0]->Integral());
0102     northx[1]->Sumw2();
0103     northx[1]->Scale(1.0/northx[1]->Integral());
0104     northx[0]->SetLineColor(kRed);
0105     northx[1]->SetLineColor(kBlue);
0106     northx[0]->SetTitle("North SMD Neutron x");
0107     northx[0]->GetYaxis()->SetTitle("Normalized Counts");
0108     northx[0]->SetStats(0);
0109     northx[0]->Draw("hist");
0110     northx[1]->Draw("same hist");
0111     TLatex* l = new TLatex;
0112     l->SetTextColor(kRed);
0113     l->DrawLatexNDC(0.7, 0.85, "Spin up");
0114     l->SetTextColor(kBlue);
0115     l->DrawLatexNDC(0.7, 0.80, "Spin down");
0116
0117     c->cd(4);
0118     southx[0]->Sumw2();
0119     southx[0]->Scale(1.0/southx[0]->Integral());
0120     southx[1]->Sumw2();
0121     southx[1]->Scale(1.0/southx[1]->Integral());
0122     southx[0]->SetLineColor(kRed);
0123     southx[1]->SetLineColor(kBlue);
0124     southx[0]->SetTitle("South SMD Neutron x");
0125     southx[0]->GetYaxis()->SetTitle("Normalized Counts");
0126     southx[0]->SetStats(0);
0127     southx[0]->Draw("hist");
0128     southx[1]->Draw("same hist");
0129     l->SetTextColor(kRed);
0130     l->DrawLatexNDC(0.7, 0.85, "Spin up");
0131     l->SetTextColor(kBlue);
0132     l->DrawLatexNDC(0.7, 0.80, "Spin down");
0133
0134     c->Update();
0135 }
0136
0137 void plot_multiplicity(TCanvas* c, TH1** hists) {
0138     c->Clear();
0139     gStyle->SetOptStat(0);
0140     c->Divide(2, 2, 0.025, 0.025);
0141
0142     c->cd(1);
0143     hists[0]->Draw();
0144     c->cd(2);
0145     hists[1]->Draw();
0146     c->cd(3);
0147     hists[2]->Draw();
0148     c->cd(4);
0149     hists[3]->Draw();
0150 }
0151
0152 void plot_asym(const char* inname, std::string outname) {
0153     std::string outname_start = outname + "(";
0154     std::string outname_end = outname + ")";
0155     TFile* f = new TFile(inname, "READ");
0156     TCanvas* c1 = new TCanvas("c1", "c1", 800, 450);
0157
0158     // hit multiplicities
0159     TH1F* h_hor_north_multiplicity = (TH1F*)f->Get("smd_hor_north_multiplicity");
0160     TH1F* h_ver_north_multiplicity = (TH1F*)f->Get("smd_ver_north_multiplicity");
0161     TH1F* h_hor_south_multiplicity = (TH1F*)f->Get("smd_hor_south_multiplicity");
0162     TH1F* h_ver_south_multiplicity = (TH1F*)f->Get("smd_ver_south_multiplicity");
0163     TH1* multiplicities[] = {h_hor_north_multiplicity, h_ver_north_multiplicity, h_hor_south_multiplicity, h_ver_south_multiplicity};
0164     plot_multiplicity(c1, multiplicities);
0165     c1->SaveAs(outname_start.c_str());
0166
0167     // beam centroid xy positions
0168     TH2F* h_xy_north = (TH2F*)f->Get("smd_xy_neutron_north");
0169     TH1F* h_ver_north_up = (TH1F*)f->Get("smd_ver_north_up");
0170     TH1F* h_ver_north_down = (TH1F*)f->Get("smd_ver_north_down");
0171     TH2F* h_xy_south = (TH2F*)f->Get("smd_xy_neutron_south");
0172     TH1F* h_ver_south_up = (TH1F*)f->Get("smd_ver_south_up");
0173     TH1F* h_ver_south_down = (TH1F*)f->Get("smd_ver_south_down");
0174     TH1* northx[] = {h_ver_north_up, h_ver_north_down};
0175     TH1* southx[] = {h_ver_south_up, h_ver_south_down};
0176     plot_xy(c1, h_xy_north, h_xy_south, northx, southx);
0177     c1->SaveAs(outname.c_str());
0178
0179     // simple asymmetry
0180     TH1F* h_phi_north_diff = (TH1F*)f->Get("smd_north_phi_diff");
0181     TH1F* h_phi_north_sum = (TH1F*)f->Get("smd_north_phi_sum");
0182     TH1F* h_phi_south_diff = (TH1F*)f->Get("smd_south_phi_diff");
0183     TH1F* h_phi_south_sum = (TH1F*)f->Get("smd_south_phi_sum");
0184
0185     // Error testing
0186     /* TH1F* h_phi_north_up = (TH1F*)f->Get("smd_north_phi_up"); */
0187     /* TH1F* h_phi_north_down = (TH1F*)f->Get("smd_north_phi_down"); */
0188     /* std::cout << "North phi up bin 1 = " << h_phi_north_up->GetBinContent(1) << ", error = " << h_phi_north_up->GetBinError(1) << std::endl; */
0189     /* std::cout << "North phi down bin 1 = " << h_phi_north_down->GetBinContent(1) << ", error = " << h_phi_north_down->GetBinError(1) << std::endl; */
0190     /* std::cout << "North phi sum bin 1 = " << h_phi_north_sum->GetBinContent(1) << ", error = " << h_phi_north_sum->GetBinError(1) << std::endl; */
0191     /* std::cout << "North phi diff bin 1 = " << h_phi_north_diff->GetBinContent(1) << ", error = " << h_phi_north_diff->GetBinError(1) << std::endl; */
0192
0193     // sqrt asymmetry
0194     TH1F* h_phi_north_L_up = (TH1F*)f->Get("smd_north_phi_L_up");
0195     TH1F* h_phi_north_L_down = (TH1F*)f->Get("smd_north_phi_L_down");
0196     TH1F* h_phi_north_R_up = (TH1F*)f->Get("smd_north_phi_R_up");
0197     TH1F* h_phi_north_R_down = (TH1F*)f->Get("smd_north_phi_R_down");
0198     TH1F* h_phi_south_L_up = (TH1F*)f->Get("smd_south_phi_L_up");
0199     TH1F* h_phi_south_L_down = (TH1F*)f->Get("smd_south_phi_L_down");
0200     TH1F* h_phi_south_R_up = (TH1F*)f->Get("smd_south_phi_R_up");
0201     TH1F* h_phi_south_R_down = (TH1F*)f->Get("smd_south_phi_R_down");
0202
0203     c1->Clear();
0204     c1->Divide(2, 2, 0.025, 0.025);
0205     TPad* p1 = (TPad*)c1->GetPad(1);
0206     TPad* p2 = (TPad*)c1->GetPad(2);
0207     TPad* p3 = (TPad*)c1->GetPad(3);
0208     TPad* p4 = (TPad*)c1->GetPad(4);
0209
0210     plot_simple(p1, h_phi_north_diff, h_phi_north_sum, "North");
0211     plot_simple(p3, h_phi_south_diff, h_phi_south_sum, "South");
0212     plot_sqrt(p2, h_phi_north_L_up, h_phi_north_R_up, h_phi_north_L_down, h_phi_north_R_down, "North");
0213     plot_sqrt(p4, h_phi_south_L_up, h_phi_south_R_up, h_phi_south_L_down, h_phi_south_R_down, "South");
0214
0215     c1->Update();
0216     c1->SaveAs(outname.c_str());
0217
0218     // Beam center correction
0219     // beam centroid xy positions
0220     h_xy_north = (TH2F*)f->Get("smd_xy_neutron_corrected_north");
0221     h_xy_south = (TH2F*)f->Get("smd_xy_neutron_corrected_south");
0222     h_ver_north_up = (TH1F*)f->Get("smd_ver_north_corrected_up");
0223     h_ver_north_down = (TH1F*)f->Get("smd_ver_north_corrected_down");
0224     h_ver_south_up = (TH1F*)f->Get("smd_ver_south_corrected_up");
0225     h_ver_south_down = (TH1F*)f->Get("smd_ver_south_corrected_down");
0226     northx[0] = h_ver_north_up;
0227     northx[1] = h_ver_north_down;
0228     southx[0] = h_ver_south_up;
0229     southx[1] = h_ver_south_down;
0230     plot_xy(c1, h_xy_north, h_xy_south, northx, southx);
0231     c1->SaveAs(outname.c_str());
0232
0233     // simple asymmetry
0234     h_phi_north_diff = (TH1F*)f->Get("smd_north_phi_diff_corrected");
0235     h_phi_north_sum = (TH1F*)f->Get("smd_north_phi_sum_corrected");
0236     h_phi_south_diff = (TH1F*)f->Get("smd_south_phi_diff_corrected");
0237     h_phi_south_sum = (TH1F*)f->Get("smd_south_phi_sum_corrected");
0238
0239     // sqrt asymmetry
0240     h_phi_north_L_up = (TH1F*)f->Get("smd_north_phi_L_up_corrected");
0241     h_phi_north_L_down = (TH1F*)f->Get("smd_north_phi_L_down_corrected");
0242     h_phi_north_R_up = (TH1F*)f->Get("smd_north_phi_R_up_corrected");
0243     h_phi_north_R_down = (TH1F*)f->Get("smd_north_phi_R_down_corrected");
0244     h_phi_south_L_up = (TH1F*)f->Get("smd_south_phi_L_up_corrected");
0245     h_phi_south_L_down = (TH1F*)f->Get("smd_south_phi_L_down_corrected");
0246     h_phi_south_R_up = (TH1F*)f->Get("smd_south_phi_R_up_corrected");
0247     h_phi_south_R_down = (TH1F*)f->Get("smd_south_phi_R_down_corrected");
0248
0249     c1->Clear();
0250     c1->Divide(2, 2, 0.025, 0.025);
0251     p1 = (TPad*)c1->GetPad(1);
0252     p2 = (TPad*)c1->GetPad(2);
0253     p3 = (TPad*)c1->GetPad(3);
0254     p4 = (TPad*)c1->GetPad(4);
0255
0256     plot_simple(p1, h_phi_north_diff, h_phi_north_sum, "Center-Corrected North");
0257     plot_simple(p3, h_phi_south_diff, h_phi_south_sum, "Center-Corrected South");
0258     plot_sqrt(p2, h_phi_north_L_up, h_phi_north_R_up, h_phi_north_L_down, h_phi_north_R_down, "Center-Corrected North");
0259     plot_sqrt(p4, h_phi_south_L_up, h_phi_south_R_up, h_phi_south_L_down, h_phi_south_R_down, "Center-Corrected South");
0260
0261     c1->Update();
0262     c1->SaveAs(outname_end.c_str());
0263
0264     return 0;
0265 }