File indexing completed on 2024-04-06 12:29:48
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031 #include <TROOT.h>
0032 #include <TChain.h>
0033 #include <TFile.h>
0034 #include <TH1D.h>
0035 #include <TH2D.h>
0036 #include <TProfile.h>
0037 #include <TProfile2D.h>
0038 #include <TFitResult.h>
0039 #include <TFitResultPtr.h>
0040 #include <TStyle.h>
0041 #include <TCanvas.h>
0042 #include <TLegend.h>
0043 #include <TPaveStats.h>
0044 #include <TPaveText.h>
0045 #include <vector>
0046 #include <string>
0047 #include <iomanip>
0048 #include <iostream>
0049 #include <fstream>
0050
0051 void makeECPlots(std::string fname = "runWithGun_Fix.root",
0052 std::string text = "Fixed Depth Calculation",
0053 std::string prefix = "Fix",
0054 bool save = false) {
0055 std::string name[4] = {"ECLL_EB", "ECLL_EBref", "ECLL_EE", "ECLL_EERef"};
0056 double xrnglo[4] = {1200.0, 1200.0, 3100.0, 3100.0};
0057 double xrnghi[4] = {1600.0, 1600.0, 3600.0, 3600.0};
0058 std::string xtitl[4] = {"R_{Cyl} (mm)", "R_{Cyl} (mm)", "|z| (mm)", "|z| (mm)"};
0059 std::string ytitl[4] = {"# X_{0} (*100)", "# X_{0} (*100)", "# X_{0} (*100)", "# X_{0} (*100)"};
0060 std::string title[4] = {"EB (Normal)", "EB (Reflected)", "EE (Normal)", "EE (Reflected)"};
0061
0062 gStyle->SetCanvasBorderMode(0);
0063 gStyle->SetCanvasColor(kWhite);
0064 gStyle->SetPadColor(kWhite);
0065 gStyle->SetFillColor(kWhite);
0066 gStyle->SetOptStat(11110);
0067 TFile *file = new TFile(fname.c_str());
0068 if (file) {
0069 char namep[100];
0070 for (int k = 0; k < 4; ++k) {
0071 TH2D *hist(0);
0072 for (int i = 0; i < 4; ++i) {
0073 if (i == 0)
0074 sprintf(namep, "%s", name[k].c_str());
0075 else
0076 sprintf(namep, "%s;%d", name[k].c_str(), i);
0077 hist = (TH2D *)file->FindObjectAny(name[k].c_str());
0078 std::cout << namep << " read out at " << hist << std::endl;
0079 if (hist != 0) {
0080 std::cout << "Entries " << hist->GetEntries() << std::endl;
0081 if (hist->GetEntries() > 0)
0082 break;
0083 }
0084 }
0085 if (hist != 0) {
0086 sprintf(namep, "%s%s", name[k].c_str(), prefix.c_str());
0087 TCanvas *pad = new TCanvas(namep, namep, 500, 500);
0088 pad->SetRightMargin(0.10);
0089 pad->SetTopMargin(0.10);
0090 hist->GetYaxis()->SetTitle(ytitl[k].c_str());
0091 hist->GetXaxis()->SetTitle(xtitl[k].c_str());
0092 hist->SetTitle(title[k].c_str());
0093 hist->GetXaxis()->SetRangeUser(xrnglo[k], xrnghi[k]);
0094 hist->GetYaxis()->SetTitleOffset(1.4);
0095 hist->Draw();
0096 pad->Update();
0097 TPaveStats *st1 = (TPaveStats *)hist->GetListOfFunctions()->FindObject("stats");
0098 if (st1 != NULL) {
0099 st1->SetY1NDC(0.70);
0100 st1->SetY2NDC(0.90);
0101 st1->SetX1NDC(0.65);
0102 st1->SetX2NDC(0.90);
0103 }
0104 TPaveText *txt1 = new TPaveText(0.50, 0.60, 0.90, 0.65, "blNDC");
0105 txt1->SetFillColor(0);
0106 txt1->AddText(text.c_str());
0107 pad->Update();
0108 if (save) {
0109 sprintf(namep, "c_%s%s.gif", name[k].c_str(), prefix.c_str());
0110 pad->Print(namep);
0111 }
0112 }
0113 }
0114 }
0115 }
0116
0117 std::vector<TCanvas *> comparePlots(std::string dirname = "EcalSimHitStudy",
0118 std::string text = "All",
0119 int mom = 10,
0120 bool ratio = false,
0121 std::string fname = "elec",
0122 bool save = false) {
0123 std::vector<TCanvas *> tcvs;
0124 gStyle->SetCanvasBorderMode(0);
0125 gStyle->SetCanvasColor(kWhite);
0126 gStyle->SetPadColor(kWhite);
0127 gStyle->SetFillColor(kWhite);
0128 gStyle->SetOptTitle(0);
0129 gStyle->SetOptFit(0);
0130 if (ratio)
0131 gStyle->SetOptStat(0);
0132 else
0133 gStyle->SetOptStat(1100);
0134
0135 std::string tags[2] = {"Old", "New"};
0136 int color[2] = {2, 4};
0137 int marker[2] = {20, 21};
0138 int style[2] = {1, 2};
0139 int rebin[16] = {50, 50, 50, 50, 2, 2, 2, 2, 2, 2, 20, 20, 20, 20, 20, 20};
0140 int type[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
0141 int edgex[16] = {0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0};
0142 std::string name1[16] = {"Etot0",
0143 "Etot1",
0144 "EtotG0",
0145 "EtotG1",
0146 "r1by250",
0147 "r1by251",
0148 "r1by90",
0149 "r1by91",
0150 "r9by250",
0151 "r9by251",
0152 "sEtaEta0",
0153 "sEtaEta1",
0154 "sEtaPhi0",
0155 "sEtaPhi1",
0156 "sPhiPhi0",
0157 "sPhiPhi1"};
0158 char name[100];
0159 TFile *file[2];
0160 TDirectory *dir[2];
0161 for (int i = 0; i < 2; ++i) {
0162 sprintf(name, "%s%d%s.root", fname.c_str(), mom, tags[i].c_str());
0163 file[i] = new TFile(name);
0164 dir[i] = (TDirectory *)file[i]->FindObjectAny(dirname.c_str());
0165 }
0166 for (int k = 0; k < 16; ++k) {
0167 TH1D *hist[2];
0168 sprintf(name, "%s", name1[k].c_str());
0169 for (int i = 0; i < 2; ++i) {
0170 hist[i] = (TH1D *)dir[i]->FindObjectAny(name);
0171 if (hist[i] != 0) {
0172 hist[i]->GetXaxis()->SetLabelOffset(0.005);
0173 hist[i]->GetXaxis()->SetTitleOffset(1.00);
0174 hist[i]->GetYaxis()->SetLabelOffset(0.005);
0175 hist[i]->GetYaxis()->SetTitleOffset(1.20);
0176 hist[i]->SetMarkerStyle(marker[i]);
0177 hist[i]->SetMarkerColor(color[i]);
0178 hist[i]->SetLineColor(color[i]);
0179 hist[i]->SetLineStyle(style[i]);
0180 hist[i]->SetLineWidth(2);
0181 }
0182 }
0183 if (hist[0] != 0 && hist[1] != 0) {
0184 double ytop(0.90), dy(0.05);
0185 double xmin = (edgex[k] == 0) ? 0.65 : 0.11;
0186 double xmin1 = (edgex[k] == 0) ? 0.55 : 0.11;
0187 double ymax = ratio ? (ytop - 0.01) : (ytop - 2 * dy - 0.01);
0188 double ymin = ratio ? (ymax - 0.045) : (ymax - 0.09);
0189 TLegend *legend = new TLegend(xmin1, ymin, xmin1 + 0.35, ymax);
0190 legend->SetFillColor(kWhite);
0191 if (ratio) {
0192 sprintf(name, "c_R%sE%d%s", name1[k].c_str(), mom, text.c_str());
0193 } else {
0194 sprintf(name, "c_%sE%d%s", name1[k].c_str(), mom, text.c_str());
0195 }
0196 TCanvas *pad = new TCanvas(name, name, 700, 500);
0197 pad->SetRightMargin(0.10);
0198 pad->SetTopMargin(0.10);
0199 if (type[k] != 0)
0200 pad->SetLogy();
0201 if (ratio) {
0202 int nbin = hist[0]->GetNbinsX();
0203 int nbinR = nbin / rebin[k];
0204 double xlow = hist[0]->GetXaxis()->GetBinLowEdge(1);
0205 double xhigh = hist[0]->GetXaxis()->GetBinLowEdge(nbin) + hist[0]->GetXaxis()->GetBinWidth(nbin);
0206 ;
0207 sprintf(name, "%sRatio", name1[k].c_str());
0208 TH1D *histr = new TH1D(name, hist[0]->GetTitle(), nbinR, xlow, xhigh);
0209 sprintf(name, "Ratio (%s/%s)", tags[0].c_str(), tags[1].c_str());
0210 histr->GetXaxis()->SetTitle(hist[0]->GetXaxis()->GetTitle());
0211 histr->GetYaxis()->SetTitle(name);
0212 histr->GetXaxis()->SetLabelOffset(0.005);
0213 histr->GetXaxis()->SetTitleOffset(1.00);
0214 histr->GetYaxis()->SetLabelOffset(0.005);
0215 histr->GetYaxis()->SetTitleOffset(1.20);
0216 histr->GetYaxis()->SetRangeUser(0.0, 2.0);
0217 histr->SetMarkerStyle(marker[0]);
0218 histr->SetMarkerColor(color[0]);
0219 histr->SetLineColor(color[0]);
0220 histr->SetLineStyle(style[0]);
0221 for (int j = 0; j < nbinR; ++j) {
0222 double tnum(0), tden(0), rnum(0), rden(0);
0223 for (int i = 0; i < rebin[k]; ++i) {
0224 int ib = j * rebin[k] + i + 1;
0225 tnum += hist[0]->GetBinContent(ib);
0226 tden += hist[1]->GetBinContent(ib);
0227 rnum += ((hist[0]->GetBinError(ib)) * (hist[0]->GetBinError(ib)));
0228 rden += ((hist[1]->GetBinError(ib)) * (hist[1]->GetBinError(ib)));
0229 }
0230 if (tden > 0 && tnum > 0) {
0231 double rat = tnum / tden;
0232 double err = rat * sqrt((rnum / (tnum * tnum)) + (rden / (tden * tden)));
0233 histr->SetBinContent(j + 1, rat);
0234 histr->SetBinError(j + 1, err);
0235 }
0236 }
0237 histr->Draw();
0238 sprintf(name, "%d GeV Electron (%s)", mom, text.c_str());
0239 legend->AddEntry(histr, name, "lp");
0240 pad->Update();
0241 TLine *line = new TLine(xlow, 1.0, xhigh, 1.0);
0242 line->SetLineColor(color[1]);
0243 line->SetLineWidth(2);
0244 line->SetLineStyle(2);
0245 line->Draw("same");
0246 pad->Update();
0247 } else {
0248 double mean[2], error[2];
0249 for (int i = 0; i < 2; ++i) {
0250 if (rebin[k] > 1)
0251 hist[i]->Rebin(rebin[k]);
0252 if (i == 0)
0253 hist[i]->Draw("hist");
0254 else
0255 hist[i]->Draw("sameshist");
0256 pad->Update();
0257 sprintf(name, "%d GeV Electron (%s %s)", mom, text.c_str(), tags[i].c_str());
0258 legend->AddEntry(hist[i], name, "lp");
0259 TPaveStats *st1 = (TPaveStats *)hist[i]->GetListOfFunctions()->FindObject("stats");
0260 if (st1 != NULL) {
0261 st1->SetLineColor(color[i]);
0262 st1->SetTextColor(color[i]);
0263 st1->SetY1NDC(ytop - dy);
0264 st1->SetY2NDC(ytop);
0265 st1->SetX1NDC(xmin);
0266 st1->SetX2NDC(xmin + 0.25);
0267 ytop -= dy;
0268 }
0269 double entries = hist[i]->GetEntries();
0270 mean[i] = hist[i]->GetMean();
0271 error[i] = (hist[i]->GetRMS()) / sqrt(entries);
0272 std::cout << text << ":" << hist[i]->GetName() << " V " << tags[i] << " Mean " << mean[i] << " RMS "
0273 << hist[i]->GetRMS() << " Error " << error[i] << std::endl;
0274 }
0275 double diff = 0.5 * (mean[0] - mean[1]) / (mean[0] + mean[1]);
0276 double ddiff =
0277 (sqrt((mean[0] * error[1]) * (mean[0] * error[1]) + (mean[1] * error[0]) * (mean[1] * error[0])) /
0278 ((mean[0] * mean[0]) + (mean[1] * mean[1])));
0279 double sign = std::abs(diff) / ddiff;
0280 std::cout << "Difference " << diff << " +- " << ddiff << " Significance " << sign << std::endl;
0281 pad->Modified();
0282 pad->Update();
0283 }
0284 if (ratio) {
0285 }
0286 legend->Draw("same");
0287 pad->Modified();
0288 pad->Update();
0289 tcvs.push_back(pad);
0290 if (save) {
0291 sprintf(name, "%s.pdf", pad->GetName());
0292 pad->Print(name);
0293 }
0294 }
0295 }
0296 return tcvs;
0297 }
