analysis_INTT/src/GeoAccepCorr.cxx

0001 #include <TCanvas.h>
0002 #include <TFile.h>
0003 #include <TH1.h>
0004 #include <TH2.h>
0005 #include <TLine.h>
0006 #include <TMath.h>
0007 #include <TObjString.h>
0008 #include <TRandom3.h>
0009 #include <TTree.h>
0010 #include <TTreeIndex.h>
0011 #include <TString.h>
0012 #include <TLegend.h>
0013 #include <TStyle.h>
0014 #include <TText.h>
0015
0016 #include <fstream>
0017 #include <iostream>
0018 #include <sstream>
0019 #include <string>
0020 #include <vector>
0021
0022 #define INCLUDE_VZ_RANGE
0023 #define INCLUDE_ETA_RANGE
0024 #include "bins.h"
0025
0026 #include "/cvmfs/sphenix.sdcc.bnl.gov/gcc-12.1.0/release/release_new/new/rootmacros/sPhenixStyle.C"
0027
0028 void convert(TH2 *h1)
0029 {
0030     TH1D *hvz = (TH1D *)h1->ProjectionY("hvz");
0031
0032     for (int i = 1; i <= h1->GetNbinsX(); ++i)
0033     {
0034         for (int j = 1; j <= h1->GetNbinsY(); ++j)
0035         {
0036             // double data_pdf = TMath::Gaus(hvz->GetBinCenter(j), -20.72, 6.390, 1);
0037             double data_pdf = TMath::Gaus(hvz->GetBinCenter(j), -4.03, 9.7, 1);
0038             if (h1->GetBinContent(i, j))
0039             {
0040                 h1->SetBinContent(i, j, data_pdf);
0041             }
0042         }
0043     }
0044
0045     delete hvz;
0046 }
0047
0048 TLine *drawline(Double_t x1, Double_t y1, Double_t x2, Double_t y2)
0049 {
0050     TLine *l = new TLine(x1, y1, x2, y2);
0051     l->SetLineColor(kRed);
0052     l->SetLineWidth(2);
0053     l->Draw();
0054     return l;
0055 }
0056
0057 void drawhoutline(TH2 *h)
0058 {
0059     auto xaxis = h->GetXaxis(), yaxis = h->GetYaxis();
0060     auto nx = xaxis->GetNbins(), ny = yaxis->GetNbins();
0061     int b[4][2] = {{0, 1}, {-1, 0}, {0, -1}, {1, 0}};
0062     for (int i = 0; i < nx; i++)
0063         for (int j = 0; j < ny; j++)
0064         {
0065             auto content = h->GetBinContent(i + 1, j + 1);
0066             for (int k = 0; k < 4; k++)
0067             {
0068                 auto neighbor = h->GetBinContent(i + 1 + b[k][0], j + 1 + b[k][1]);
0069                 if ((content || neighbor) && !(content && neighbor))
0070                 {
0071                     auto cx = xaxis->GetBinCenter(i + 1), wx = xaxis->GetBinWidth(i + 1), cy = yaxis->GetBinCenter(j + 1), wy = yaxis->GetBinWidth(j + 1), fx = (double)b[k][0], fy = (double)b[k][1];
0072                     drawline(cx + wx * (fx + fy) / 2., cy + wy * (fy + fx) / 2, cx + wx * (fx - fy) / 2., cy + wy * (fy - fx) / 2.);
0073                 }
0074             }
0075         }
0076 }
0077
0078 int main(int argc, char *argv[])
0079 {
0080     if (argc != 4)
0081     {
0082         std::cout << "Usage: ./GeoAccepCorr [datatreefile (minitree)] [simtreefile (minitree)] [outfilename]" << std::endl;
0083         return 0;
0084     }
0085
0086     for (int i = 0; i < argc; i++)
0087     {
0088         std::cout << "argv[" << i << "] = " << argv[i] << std::endl;
0089     }
0090
0091     SetsPhenixStyle();
0092
0093     TString datatreefile = TString(argv[1]);
0094     TString simtreefile = TString(argv[2]);
0095     TString outfilename = TString(argv[3]);
0096
0097     TFile *fdata = new TFile(datatreefile, "READ");
0098     TTree *tdata = (TTree *)fdata->Get("minitree");
0099
0100     TFile *fsim = new TFile(simtreefile, "READ");
0101     TTree *tsim = (TTree *)fsim->Get("minitree");
0102
0103     int nfeta = neta * 100;
0104     int nfvz = nvz * 100;
0105
0106     auto evtsel = "PV_z>=-10&&PV_z<=10";
0107
0108     TH2D *hM_data = new TH2D("hM_data", "hM_data", nfeta, etamin, etamax, nfvz, vzmin, vzmax);
0109     tdata->Project("hM_data", "PV_z:recotklraw_eta", evtsel, "");
0110     convert(hM_data);
0111     TH2D *hM_data_coarse = (TH2D *)hM_data->Clone("hM_data_coarse");
0112     hM_data_coarse->RebinX(nfeta / neta);
0113     hM_data_coarse->RebinY(nfvz / nvz);
0114
0115     TH2D *hM_sim = new TH2D("hM_sim", "hM_sim", nfeta, etamin, etamax, nfvz, vzmin, vzmax);
0116     tsim->Project("hM_sim", "PV_z:recotklraw_eta", evtsel, "");
0117     convert(hM_sim);
0118     TH2D *hM_sim_coarse = (TH2D *)hM_sim->Clone("hM_sim_coarse");
0119     hM_sim_coarse->RebinX(nfeta / neta);
0120     hM_sim_coarse->RebinY(nfvz / nvz);
0121
0122     TH2D *hM_ratio = (TH2D *)hM_sim_coarse->Clone("hM_ratio");
0123     hM_ratio->Divide(hM_data_coarse);
0124     // hM_ratio->SetStats(0);
0125     // TH2D *hM_ratio = new TH2D("hM_ratio", "hM_ratio", neta, etamin, etamax, nvz, vzmin, vzmax);
0126
0127     auto ext_accep_map = (TH2D *)hM_ratio->Clone("ext_accep_map");
0128     for (int i = 1; i <= neta; i++)
0129     {
0130         for (int j = 1; j <= nvz; j++)
0131         {
0132             // calculate uncertainty
0133             double data = hM_data_coarse->GetBinContent(i, j);
0134             double sim = hM_sim_coarse->GetBinContent(i, j);
0135             double data_err = hM_data_coarse->GetBinError(i, j);
0136             double sim_err = hM_sim_coarse->GetBinError(i, j);
0137             double ratio = -1E9;
0138             double ratio_err = -1E9;
0139             if (data == 0 || sim == 0)
0140             {
0141                 ratio = 0;
0142                 ratio_err = 0;
0143             }
0144             else
0145             {
0146                 ratio = sim / data;
0147                 ratio_err = ratio * sqrt(pow(sim_err / sim, 2) + pow(data_err / data, 2));
0148             }
0149             // hM_ratio->SetBinContent(i, j, ratio);
0150             hM_ratio->SetBinError(i, j, ratio_err);
0151             std::cout << "Bin (eta, vtxz) = (" << hM_data_coarse->GetXaxis()->GetBinCenter(i) << ", " << hM_data_coarse->GetYaxis()->GetBinCenter(j) << "): data = " << data << " +/- " << data_err << ", sim = " << sim << " +/- " << sim_err << ", ratio = " << ratio << " +/- " << ratio_err << std::endl;
0152
0153             if (hM_ratio->GetBinContent(i, j) < 0.75 || hM_ratio->GetBinContent(i, j) > 1.25)
0154             {
0155                 ext_accep_map->SetBinContent(i, j, 0);
0156             }
0157             else
0158             {
0159                 ext_accep_map->SetBinContent(i, j, 1);
0160             }
0161             ext_accep_map->SetBinError(i, j, 0);
0162         }
0163     }
0164
0165     // print out for acceptance.h
0166     for (int z = nvz; z >= 1; z--)
0167     {
0168         for (int x = 1; x <= neta; x++)
0169         {
0170             std::cout << ext_accep_map->GetBinContent(x, z) << ",";
0171         }
0172         std::cout << std::endl;
0173     }
0174
0175     TCanvas *c = new TCanvas("c", "c", 800, 700);
0176     gPad->SetRightMargin(0.16);
0177     gPad->SetTopMargin(0.08);
0178     c->cd();
0179     hM_ratio->GetXaxis()->SetTitle("#eta");
0180     hM_ratio->GetYaxis()->SetTitle("vtx_{Z} [cm]");
0181     hM_ratio->GetZaxis()->SetRangeUser(hM_ratio->GetMinimum(0)*0.9, hM_ratio->GetMaximum()*1.1);
0182     gStyle->SetPaintTextFormat("1.3f");
0183     hM_ratio->SetContour(1000);
0184     hM_ratio->SetMarkerSize(0.5);
0185     hM_ratio->Draw("colztext");
0186     drawhoutline(ext_accep_map);
0187     c->SaveAs(outfilename+".pdf");
0188
0189     c->Clear();
0190     c->cd();
0191     hM_data->GetXaxis()->SetTitle("#eta");
0192     hM_data->GetYaxis()->SetTitle("vtx_{Z} [cm]");
0193     hM_data->SetContour(1000);
0194     hM_data->Draw("colz");
0195     // TLegend *l = new TLegend(0.2, 0.2, 0.4, 0.4);
0196     // l->SetTextSize(0.04);
0197     // l->SetFillStyle(0);
0198     // l->AddEntry("", "Data", "");
0199     // l->Draw();
0200     TText *t = new TText();
0201     // right and bottom adjusted
0202     t->SetTextAlign(31);
0203     t->SetTextSize(0.04);
0204     t->DrawTextNDC(1-gPad->GetRightMargin(), (1 - gPad->GetTopMargin()) + 0.01, "Data");
0205     c->SaveAs(outfilename+"_hM_data.pdf");
0206     c->SaveAs(outfilename+"_hM_data.png");
0207
0208     // l->Clear();
0209     c->Clear();
0210     c->cd();
0211     hM_sim->GetXaxis()->SetTitle("#eta");
0212     hM_sim->GetYaxis()->SetTitle("vtx_{Z} [cm]");
0213     hM_sim->SetContour(1000);
0214     hM_sim->Draw("colz");
0215     // l->AddEntry("", "Simulation", "");
0216     // l->Draw();
0217     t->DrawTextNDC(1-gPad->GetRightMargin(), (1 - gPad->GetTopMargin()) + 0.01, "Simulation");
0218     c->SaveAs(outfilename+"_hM_sim.pdf");
0219     c->SaveAs(outfilename+"_hM_sim.png");
0220
0221     // l->Clear();
0222     c->Clear();
0223     c->cd();
0224     hM_data_coarse->GetXaxis()->SetTitle("#eta");
0225     hM_data_coarse->GetYaxis()->SetTitle("vtx_{Z} [cm]");
0226     hM_data_coarse->SetContour(1000);
0227     hM_data_coarse->Draw("colz");
0228     // l->AddEntry("", "Data", "");
0229     // l->Draw();
0230     t->DrawTextNDC(1-gPad->GetRightMargin(), (1 - gPad->GetTopMargin()) + 0.01, "Data");
0231     c->SaveAs(outfilename+"_hM_data_coarse.pdf");
0232
0233     // l->Clear();
0234     c->Clear();
0235     c->cd();
0236     hM_sim_coarse->GetXaxis()->SetTitle("#eta");
0237     hM_sim_coarse->GetYaxis()->SetTitle("vtx_{Z} [cm]");
0238     hM_sim_coarse->SetContour(1000);
0239     hM_sim_coarse->Draw("colz");
0240     // l->AddEntry("", "Simulation", "");
0241     // l->Draw();
0242     t->DrawTextNDC(1-gPad->GetRightMargin(), (1 - gPad->GetTopMargin()) + 0.01, "Simulation");
0243     c->SaveAs(outfilename+"_hM_sim_coarse.pdf");
0244
0245     TFile *fout = new TFile(outfilename+".root", "RECREATE");
0246     fout->cd();
0247     c->Write();
0248     hM_data->Write();
0249     hM_sim->Write();
0250     hM_data_coarse->Write();
0251     hM_sim_coarse->Write();
0252     hM_ratio->Write();
0253     ext_accep_map->Write();
0254     fout->Close();
0255 }