Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:28:59

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