Calorimetry/test/noise_plots.C

0001 void noise_plots(char* newFile="DQM_fastsim.root",char* refFile="DQM_fullsim.root") {
0002
0003   gROOT ->Reset();
0004   gROOT ->SetBatch();
0005
0006   //=========  settings ====================
0007   gROOT->SetStyle("Plain");
0008   gStyle->SetPadGridX(kTRUE);
0009   gStyle->SetPadGridY(kTRUE);
0010   gStyle->SetPadRightMargin(0.07);
0011   gStyle->SetPadLeftMargin(0.13);
0012   char* refLabel("reference histogram");
0013   char* newLabel("new histogram");
0014
0015   delete gROOT->GetListOfFiles()->FindObject(refFile);
0016   delete gROOT->GetListOfFiles()->FindObject(newFile);
0017
0018   TText* te = new TText();
0019   TFile * sfile = new TFile(newFile);
0020   TDirectory * sdir=gDirectory;
0021   TFile * rfile = new TFile(refFile);
0022   TDirectory * rdir=gDirectory;
0023
0024   if(sfile->GetDirectory("DQMData/Run 1/HcalRecHitsV")) sfile->cd("DQMData/Run 1/HcalRecHitsV/Run summary/HcalRecHitTask");
0025   else if(sfile->GetDirectory("DQMData/HcalRecHitsV/HcalRecHitTask"))sfile->cd("DQMData/HcalRecHitsV/HcalRecHitTask");
0026   sdir=gDirectory;
0027   TList *sl= sdir->GetListOfKeys();
0028
0029   if(rfile->GetDirectory("DQMData/Run 1/HcalRecHitsV")) rfile->cd("DQMData/Run 1/HcalRecHitsV/Run summary/HcalRecHitTask");
0030   else if(rfile->GetDirectory("DQMData/HcalRecHitsV/HcalRecHitTask"))rfile->cd("DQMData/HcalRecHitsV/HcalRecHitTask");
0031   rdir=gDirectory;
0032   TList *rl= rdir->GetListOfKeys();
0033
0034   TCanvas *canvas;
0035
0036   TH1F *sh1,*rh1;
0037   TH1F *sh2,*rh2;
0038   TH1F *sh3,*rh3;
0039   TH1F *sh4,*rh4;
0040   TH1F *sh5,*rh5;
0041
0042
0043   rdir->GetObject("e_hb",rh1);
0044   sdir->GetObject("e_hb",sh1);
0045   rdir->GetObject("e_he",rh2);
0046   sdir->GetObject("e_he",sh2);
0047   rdir->GetObject("e_hfl",rh3);
0048   sdir->GetObject("e_hfl",sh3);
0049   rdir->GetObject("e_hfs",rh4);
0050   sdir->GetObject("e_hfs",sh4);
0051   rdir->GetObject("e_ho",rh5);
0052   sdir->GetObject("e_ho",sh5);
0053
0054   rh1->GetXaxis()->SetTitle("E");
0055   rh1->GetYaxis()->SetTitleSize(0.05);
0056   rh1->GetYaxis()->SetTitleOffset(1.2);
0057   //  rh1->GetYaxis()->SetRangeUser(0.4,1.);
0058   //  sh1->GetYaxis()->SetRangeUser(0.4,1.);
0059   //  rh1->Rebin(5);
0060   //  sh1->Rebin(5);
0061   NormalizeHistograms(rh1,sh1);
0062   double maxH = max(rh1,sh1)*1.1;
0063   rh1->GetYaxis()->SetRangeUser(0.,maxH);
0064   sh1->GetYaxis()->SetRangeUser(0.,maxH);
0065
0066   rh2->GetXaxis()->SetTitle("E");
0067   rh2->GetYaxis()->SetTitleSize(0.05);
0068   rh2->GetYaxis()->SetTitleOffset(1.2);
0069   //  rh2->GetYaxis()->SetRangeUser(0.4,1.);
0070   //  sh2->GetYaxis()->SetRangeUser(0.4,1.);
0071   //  rh2->Rebin(5);
0072   //  sh2->Rebin(5);
0073   NormalizeHistograms(rh2,sh2);
0074   double maxH = max(rh2,sh2)*1.1;
0075   rh2->GetYaxis()->SetRangeUser(0.,maxH);
0076   sh2->GetYaxis()->SetRangeUser(0.,maxH);
0077
0078   rh3->GetXaxis()->SetTitle("E");
0079   rh3->GetYaxis()->SetTitleSize(0.05);
0080   rh3->GetYaxis()->SetTitleOffset(1.2);
0081   NormalizeHistograms(rh3,sh3);
0082   double maxH = max(rh3,sh3)*1.1;
0083   rh3->GetYaxis()->SetRangeUser(0.,maxH);
0084   sh3->GetYaxis()->SetRangeUser(0.,maxH);
0085
0086   rh4->GetXaxis()->SetTitle("E");
0087   rh4->GetYaxis()->SetTitleSize(0.05);
0088   rh4->GetYaxis()->SetTitleOffset(1.2);
0089   NormalizeHistograms(rh4,sh4);
0090   double maxH = max(rh4,sh4)*1.1;
0091   rh4->GetYaxis()->SetRangeUser(0.,maxH);
0092   sh4->GetYaxis()->SetRangeUser(0.,maxH);
0093
0094   rh5->GetXaxis()->SetTitle("E");
0095   rh5->GetYaxis()->SetTitleSize(0.05);
0096   rh5->GetYaxis()->SetTitleOffset(1.2);
0097   NormalizeHistograms(rh5,sh5);
0098   double maxH = max(rh5,sh5)*1.1;
0099   rh5->GetYaxis()->SetRangeUser(0.,maxH);
0100   sh5->GetYaxis()->SetRangeUser(0.,maxH);
0101
0102   canvas = new TCanvas("Noise","Noise",1000,1400);
0103
0104
0105   TH1F * r[5]={rh1,rh2,rh3,rh4,rh5};
0106   TH1F * s[5]={sh1,sh2,sh3,sh4,sh5};
0107   plotBuilding(canvas,s, r,5,
0108            te,"UU",-1);//, 1, true, 0);
0109
0110   canvas->cd();
0111   l = new TLegend(0.50,0.14,0.90,0.19);
0112   l->SetTextSize(0.016);
0113   l->SetLineColor(1);
0114   l->SetLineWidth(1);
0115   l->SetLineStyle(1);
0116   l->SetFillColor(0);
0117   l->SetBorderSize(3);
0118   l->AddEntry(rh1,refLabel,"LPF");
0119   l->AddEntry(sh1,newLabel,"LPF");
0120   l->Draw();
0121   TString namepdf = "hcal_noise.pdf";
0122   canvas->Print(namepdf);
0123   delete l;
0124
0125   delete canvas;
0126
0127 }
0128
0129 void plotBuilding(TCanvas *canvas, TH1F **s, TH1F **r, int n,TText* te,
0130           char * option, double startingY, double startingX = .1,bool fit = false, unsigned int logx=0){
0131   canvas->Divide(2,(n+1)/2); //this should work also for odd n
0132   for(int i=0; i<n;i++){
0133     s[i]->SetMarkerStyle(20);
0134     r[i]->SetMarkerStyle(21);
0135     s[i]->SetMarkerColor(2);
0136     r[i]->SetMarkerColor(4);
0137     s[i]->SetMarkerSize(0.7);
0138     r[i]->SetMarkerSize(0.7);
0139     s[i]->SetLineColor(1);
0140     r[i]->SetLineColor(1);
0141     s[i]->SetLineWidth(1);
0142     r[i]->SetLineWidth(1);
0143
0144     TPad *pad=canvas->cd(i+1);
0145     setStats(s[i],r[i], -1, 0, fit);
0146     if((logx>>i)&1)pad->SetLogx();
0147     r[i]->Draw("lpf");
0148     s[i]->Draw("sameslpf");
0149   }
0150
0151 }
0152
0153 void setStats(TH1* s,TH1* r, double startingY, double startingX = .1,bool fit){
0154   if (startingY<0){
0155     s->SetStats(0);
0156     r->SetStats(0);
0157   } else {
0158     //gStyle->SetOptStat(1001);
0159     s->SetStats(1);
0160     r->SetStats(1);
0161
0162     if (fit){
0163       s->Fit("gaus");
0164       TF1* f1 = (TF1*) s->GetListOfFunctions()->FindObject("gaus");
0165       f1->SetLineColor(2);
0166       f1->SetLineWidth(1);
0167     }
0168     s->Draw();
0169     gPad->Update();
0170     TPaveStats* st1 = (TPaveStats*) s->GetListOfFunctions()->FindObject("stats");
0171     if (fit) {st1->SetOptFit(0010);    st1->SetOptStat(1001);}
0172     st1->SetX1NDC(startingX);
0173     st1->SetX2NDC(startingX+0.30);
0174     st1->SetY1NDC(startingY+0.20);
0175     st1->SetY2NDC(startingY+0.35);
0176     st1->SetTextColor(2);
0177     if (fit) {
0178       r->Fit("gaus");
0179       TF1* f2 = (TF1*) r->GetListOfFunctions()->FindObject("gaus");
0180       f2->SetLineColor(4);
0181       f2->SetLineWidth(1);
0182     }
0183     r->Draw();
0184     gPad->Update();
0185     TPaveStats* st2 = (TPaveStats*) r->GetListOfFunctions()->FindObject("stats");
0186     if (fit) {st2->SetOptFit(0010);    st2->SetOptStat(1001);}
0187     st2->SetX1NDC(startingX);
0188     st2->SetX2NDC(startingX+0.30);
0189     st2->SetY1NDC(startingY);
0190     st2->SetY2NDC(startingY+0.15);
0191     st2->SetTextColor(4);
0192   }
0193 }
0194
0195 void NormalizeHistograms(TH1F* h1, TH1F* h2)
0196 {
0197   if (h1==0 || h2==0) return;
0198   float scale1 = -9999.9;
0199   float scale2 = -9999.9;
0200
0201   if ( h1->Integral() != 0 && h2->Integral() != 0 ){
0202     scale1 = 1.0/(float)h1->Integral();
0203     scale2 = 1.0/(float)h2->Integral();
0204
0205     //h1->Sumw2();
0206     //h2->Sumw2();
0207     h1->Scale(scale1);
0208     h2->Scale(scale2);
0209
0210   }
0211 }
0212
0213 double max(TH1F* h1, TH1F* h2)
0214 {
0215   double a = h1->GetMaximum();
0216   double b = h2->GetMaximum();
0217   return (a > b) ? a : b;
0218 }