File indexing completed on 2023-03-17 10:38:31
0001 #include "DrawPlot.h"
0002
0003
0004 #include <iostream>
0005 #include <sstream>
0006 #include <iomanip>
0007
0008 #include <cmath>
0009
0010 #include "TDirectory.h"
0011 #include "TCanvas.h"
0012 #include "TLine.h"
0013 #include "TBranch.h"
0014 #include "TPaveStats.h"
0015 #include "TStyle.h"
0016 #include "TGaxis.h"
0017
0018
0019
0020
0021 DrawPlot::DrawPlot(const unsigned int iterationNumber, const bool summaryFile):
0022 outpath_(nullptr), file_(nullptr), fileZeroApe_(nullptr), designFile_(nullptr), baselineTreeX_(nullptr), baselineTreeY_(nullptr), delta0_(nullptr),
0023 legendEntry_("data (final APE)"), legendEntryZeroApe_("data (APE=0)"), designLegendEntry_("MCideal"),
0024 legendXmin_(0.41), legendYmin_(0.27), legendXmax_(0.71), legendYmax_(0.42),
0025 thesisMode_(false)
0026 {
0027 std::stringstream ss_inpath, ss_inpathZeroApe;
0028 ss_inpath<<"$CMSSW_BASE/src/Alignment/APEEstimation/hists/workingArea/iter";
0029 ss_inpathZeroApe<<ss_inpath.str()<<"0/";
0030 ss_inpath<<iterationNumber<<"/";
0031
0032 const TString* inpath = new TString(ss_inpath.str().c_str());
0033 const TString* inpathZeroApe = new TString(ss_inpathZeroApe.str().c_str());
0034 outpath_ = new TString(inpath->Copy().Append("plots/"));
0035 const TString rootFileName(summaryFile ? "allData_resultsFile.root" : "allData.root");
0036 const TString* fileName = new TString(inpath->Copy().Append(rootFileName));
0037 const TString* fileNameZeroApe = new TString(inpathZeroApe->Copy().Append(rootFileName));
0038 const TString* designFileName = new TString("$CMSSW_BASE/src/Alignment/APEEstimation/hists/Design/baseline/" + rootFileName);
0039 const TString* baselineFileName = new TString("$CMSSW_BASE/src/Alignment/APEEstimation/hists/Design/baseline/allData_baselineApe.root");
0040
0041 std::cout<<"\n";
0042 std::cout<<"Outpath: "<<*outpath_<<"\n";
0043 std::cout<<"File name (final APE): "<<*fileName<<"\n";
0044 std::cout<<"File name (zero APE): "<<*fileNameZeroApe<<"\n";
0045 std::cout<<"Design file name: "<<*designFileName<<"\n";
0046 std::cout<<"Baseline file name: "<<*baselineFileName<<"\n";
0047 std::cout<<"\n";
0048
0049 if(iterationNumber!=0)file_ = new TFile(*fileName, "READ");
0050 fileZeroApe_ = new TFile(*fileNameZeroApe, "READ");
0051 designFile_ = new TFile(*designFileName, "READ");
0052 TFile* baselineFile = new TFile(*baselineFileName, "READ");
0053
0054
0055
0056
0057
0058
0059 if(baselineFile){
0060 baselineFile->GetObject("iterTreeX", baselineTreeX_);
0061 baselineFile->GetObject("iterTreeY", baselineTreeY_);
0062 }
0063
0064 if(!baselineTreeX_)std::cout<<"Baseline tree for x coordinate not found, cannot draw baselines!\n";
0065 if(!baselineTreeY_)std::cout<<"Baseline tree for y coordinate not found, cannot draw baselines!\n";
0066
0067 baselineTreeX_->SetDirectory(nullptr);
0068 baselineTreeY_->SetDirectory(nullptr);
0069
0070 delete inpath;
0071 delete fileName;
0072 delete fileNameZeroApe;
0073 delete designFileName;
0074 delete baselineFileName;
0075 }
0076
0077
0078
0079 DrawPlot::~DrawPlot(){
0080 if(outpath_)delete outpath_;
0081 if(file_)file_->Close();
0082 if(fileZeroApe_)fileZeroApe_->Close();
0083 if(designFile_)designFile_->Close();
0084 if(baselineTreeX_)baselineTreeX_->Delete();
0085 if(baselineTreeY_)baselineTreeY_->Delete();
0086 }
0087
0088
0089
0090 void
0091 DrawPlot::setLegendEntry(const TString& legendEntry, const TString& legendEntryZeroApe, const TString& designLegendEntry){
0092 legendEntry_ = legendEntry;
0093 legendEntryZeroApe_ = legendEntryZeroApe;
0094 designLegendEntry_ = designLegendEntry;
0095
0096 if(thesisMode_){
0097
0098 }
0099 }
0100
0101
0102 void
0103 DrawPlot::setLegendCoordinate(const double legendXmin, const double legendYmin, const double legendXmax, const double legendYmax){
0104 legendXmin_ = legendXmin;
0105 legendYmin_ = legendYmin;
0106 legendXmax_ = legendXmax;
0107 legendYmax_ = legendYmax;
0108 }
0109
0110
0111 DrawPlot::LegendEntries
0112 DrawPlot::adjustLegendEntry(const TString& histName, TH1*& hist, TH1*& histZeroApe, TH1*& designHist){
0113 LegendEntries legendEntries;
0114 legendEntries.legendEntry = legendEntry_;
0115 legendEntries.legendEntryZeroApe = legendEntryZeroApe_;
0116 legendEntries.designLegendEntry = designLegendEntry_;
0117 if(!thesisMode_)return legendEntries;
0118
0119 double mean(-999.);
0120 double meanZeroApe(-999.);
0121 double meanDesign(-999.);
0122 double rms(-999.);
0123 double rmsZeroApe(-999.);
0124 double rmsDesign(-999.);
0125 if(hist){
0126 mean = hist->GetMean();
0127 rms = hist->GetRMS();
0128 hist->SetTitle("");
0129 }
0130 if(histZeroApe){
0131 meanZeroApe = histZeroApe->GetMean();
0132 rmsZeroApe = histZeroApe->GetRMS();
0133 histZeroApe->SetTitle("");
0134 }
0135 if(designHist){
0136 meanDesign = designHist->GetMean();
0137 rmsDesign = designHist->GetRMS();
0138 designHist->SetTitle("");
0139 }
0140
0141
0142 std::string mode("");
0143 unsigned int precision(0);
0144 std::string unit("");
0145
0146 if(histName.Contains("h_norChi2")){
0147 mode = "mean";
0148 precision = 2;
0149 unit = "";
0150 }
0151 else if(histName.Contains("h_etaSig")){
0152 mode = "rms";
0153 precision = 0;
0154 unit = "";
0155 }
0156 else if(histName.Contains("h_etaErr")){
0157 mode = "mean";
0158 precision = 2;
0159 unit = " x10^{-3}";
0160
0161 mean *= 1000.;
0162 meanZeroApe *= 1000.;
0163 meanDesign *= 1000.;
0164 }
0165 else if(histName.Contains("h_phiSig")){
0166 mode = "rms";
0167 precision = 0;
0168 unit = "";
0169 }
0170 else if(histName.Contains("h_phiErr")){
0171 mode = "mean";
0172 precision = 1;
0173 unit = " x10^{-3} ^{o}";
0174
0175 mean *= 1000.;
0176 meanZeroApe *= 1000.;
0177 meanDesign *= 1000.;
0178 }
0179 else if(histName.Contains("h_phi")){
0180
0181
0182
0183 }
0184 else if(histName.Contains("h_ptSig")){
0185 mode = "rms";
0186 precision = 1;
0187 unit = "";
0188 }
0189 else if(histName.Contains("h_ptErr")){
0190 mode = "mean";
0191 precision = 2;
0192 unit = " GeV";
0193 }
0194 else if(histName.Contains("h_pt")){
0195 mode = "mean";
0196 precision = 1;
0197 unit = " GeV";
0198 }
0199 else if(histName.Contains("h_prob")){
0200 mode = "mean";
0201 precision = 2;
0202 unit = "";
0203 }
0204 else if(histName.Contains("h_p")){
0205 mode = "mean";
0206 precision = 1;
0207 unit = " GeV";
0208 }
0209 else if(histName.Contains("h_d0BeamspotSig")){
0210 mode = "rms";
0211 precision = 3;
0212 unit = "";
0213 }
0214 else if(histName.Contains("h_d0BeamspotErr")){
0215 mode = "mean";
0216 precision = 1;
0217 unit = " #mum";
0218
0219 mean *= 10000.;
0220 meanZeroApe *= 10000.;
0221 meanDesign *= 10000.;
0222 }
0223 else if(histName.Contains("h_d0Beamspot")){
0224 mode = "rms";
0225 precision = 1;
0226 unit = " #mum";
0227
0228 rms *= 10000.;
0229 rmsZeroApe *= 10000.;
0230 rmsDesign *= 10000.;
0231 }
0232 else if(histName.Contains("h_dzSig")){
0233 mode = "rms";
0234 precision = 0;
0235 unit = "";
0236 }
0237 else if(histName.Contains("h_dzErr")){
0238 mode = "mean";
0239 precision = 1;
0240 unit = " #mum";
0241
0242 mean *= 10000.;
0243 meanZeroApe *= 10000.;
0244 meanDesign *= 10000.;
0245 }
0246 else if(histName.Contains("h_dz")){
0247 mode = "rms";
0248 precision = 1;
0249 unit = " cm";
0250 }
0251 else if(histName.Contains("h_NorResX") || histName.Contains("h_NorResY")){
0252 mode = "rms";
0253 precision = 2;
0254 unit = "";
0255 }
0256 else if(histName.Contains("h_ResX") || histName.Contains("h_ResY")){
0257 mode = "rms";
0258 precision = 1;
0259 unit = " #mum";
0260 }
0261
0262 std::stringstream legendEntry;
0263 std::stringstream legendEntryZeroApe;
0264 std::stringstream designLegendEntry;
0265 if(mode == "mean"){
0266 legendEntry<<" #mu="<<std::fixed<<std::setprecision(precision)<<mean<<unit;
0267 legendEntryZeroApe<<" #mu="<<std::fixed<<std::setprecision(precision)<<meanZeroApe<<unit;
0268 designLegendEntry<<" #mu="<<std::fixed<<std::setprecision(precision)<<meanDesign<<unit;
0269 }
0270 else if(mode == "rms"){
0271 legendEntry<<" rms="<<std::fixed<<std::setprecision(precision)<<rms<<unit;
0272 legendEntryZeroApe<<" rms="<<std::fixed<<std::setprecision(precision)<<rmsZeroApe<<unit;
0273 designLegendEntry<<" rms="<<std::fixed<<std::setprecision(precision)<<rmsDesign<<unit;
0274 }
0275 else if(mode == ""){;}
0276 else{
0277 std::cout<<"Incorrect mode for legend adjustment, skip adjustment: "<<mode<<"\n";
0278 return legendEntries;
0279 }
0280
0281 legendEntries.legendEntry += legendEntry.str().c_str();
0282 legendEntries.legendEntryZeroApe += legendEntryZeroApe.str().c_str();
0283 legendEntries.designLegendEntry += designLegendEntry.str().c_str();
0284
0285 return legendEntries;
0286 }
0287
0288
0289 void
0290 DrawPlot::drawPlot(const TString& pluginName, const TString& histName, const bool normalise, const bool plotZeroApe){
0291 TString* plugin = new TString(pluginName.Copy());
0292 if(!plugin->IsNull())plugin->Append("/");
0293 for(unsigned int iSector=1; ; ++iSector){
0294 std::stringstream ss_sectorName, ss_sector;
0295 ss_sectorName<<"Sector_"<<iSector;
0296 ss_sector<<*plugin<<ss_sectorName.str()<<"/";
0297 TDirectory* dir(nullptr);
0298
0299 dir = (TDirectory*)designFile_->TDirectory::GetDirectory(ss_sector.str().c_str());
0300 if(!dir)break;
0301
0302 TH1* SectorName(nullptr);
0303 designFile_->GetObject((ss_sector.str()+"z_name;1").c_str(), SectorName);
0304 const TString sectorName(SectorName ? SectorName->GetTitle() : ss_sectorName.str().c_str());
0305
0306
0307 TTree* baselineTree(nullptr);
0308 if(histName=="h_residualWidthX1"){baselineTree = baselineTreeX_;}
0309 else if(histName=="h_residualWidthY1"){baselineTree = baselineTreeY_;}
0310 if(baselineTree){
0311 std::stringstream ss_branch;
0312 ss_branch<<"Ape_Sector_"<<iSector;
0313 TBranch* branch(nullptr);
0314 branch = baselineTree->GetBranch(ss_branch.str().c_str());
0315 if(branch){
0316 double delta0(999.);
0317 branch->SetAddress(&delta0);
0318 branch->GetEntry(0);
0319 delta0_ = new double(std::sqrt(delta0));
0320 }
0321 else delta0_ = 0;
0322 }
0323 else delta0_ = 0;
0324
0325 if(histName=="h_entriesX" || histName=="h_entriesY" ||
0326 histName=="h_ResX" || histName=="h_ResY" ||
0327 histName=="h_NorResX" || histName=="h_NorResY")ss_sector<<"Results/";
0328
0329 ss_sector<<histName;
0330 const TString fullName(ss_sector.str().c_str());
0331 this->printHist(fullName, histName.Copy().Append("_").Append(sectorName), normalise, plotZeroApe);
0332
0333 if(delta0_)delete delta0_;
0334 }
0335 }
0336
0337
0338
0339 void
0340 DrawPlot::drawTrackPlot(const TString& pluginName, const TString& histName, const bool normalise, const bool plotZeroApe){
0341 TString* plugin = new TString(pluginName.Copy());
0342 if(!plugin->IsNull())plugin->Append("/");
0343 std::stringstream ss_sector;
0344 ss_sector<<*plugin<<"TrackVariables/"<<histName;
0345 const TString fullName(ss_sector.str().c_str());
0346 this->printHist(fullName, histName, normalise, plotZeroApe);
0347 }
0348
0349
0350
0351 void
0352 DrawPlot::drawEventPlot(const TString& pluginName, const TString& histName, const bool normalise, const bool plotZeroApe){
0353 TString* plugin = new TString(pluginName.Copy());
0354 if(!plugin->IsNull())plugin->Append("/");
0355 std::stringstream ss_sector;
0356 ss_sector<<*plugin<<"EventVariables/"<<histName;
0357 const TString fullName(ss_sector.str().c_str());
0358 this->printHist(fullName, histName, normalise, plotZeroApe);
0359 }
0360
0361
0362
0363 void
0364 DrawPlot::printHist(const TString& fullName, const TString& sectorName, const bool normalise, const bool plotZeroApe){
0365 if(thesisMode_)gStyle->SetOptStat(0);
0366
0367 TH1* hist(nullptr);
0368 TH1* histZeroApe(nullptr);
0369 TH1* designHist(nullptr);
0370 if(file_)file_->GetObject(fullName + ";1", hist);
0371 if(fileZeroApe_)fileZeroApe_->GetObject(fullName + ";1", histZeroApe);
0372 designFile_->GetObject(fullName + ";1", designHist);
0373 if(hist && !plotZeroApe)histZeroApe = 0;
0374 if(!(hist || histZeroApe) || !designHist){std::cout<<"Histogram not found in file: "<<fullName<<"\n"; return;}
0375 else std::cout<<"Drawing histogram: "<<fullName<<"\n";
0376
0377 std::vector<TH1*> v_hist;
0378 v_hist.push_back(designHist);
0379 if(histZeroApe)v_hist.push_back(histZeroApe);
0380 if(hist)v_hist.push_back(hist);
0381
0382 if(normalise)this->scale(v_hist, 100.);
0383
0384 const double maxY(this->maximumY(v_hist));
0385
0386 this->setRangeUser(v_hist, 0., 1.1*maxY);
0387
0388
0389 const double maxYAxis(1.1*maxY);
0390 if(maxYAxis<10.)TGaxis::SetMaxDigits(3);
0391 else TGaxis::SetMaxDigits(4);
0392 if(sectorName.Contains("h_weightX") || sectorName.Contains("h_weightY")){
0393 if(maxYAxis<0.6){
0394 if(hist)hist->SetNdivisions(506, "Y");
0395 if(histZeroApe)histZeroApe->SetNdivisions(506, "Y");
0396 if(designHist)designHist->SetNdivisions(506, "Y");
0397 }
0398 if(maxYAxis<0.3){
0399 if(hist)hist->SetNdivisions(503, "Y");
0400 if(histZeroApe)histZeroApe->SetNdivisions(503, "Y");
0401 if(designHist)designHist->SetNdivisions(503, "Y");
0402 }
0403 }
0404 if(sectorName.Contains("h_d0BeamspotErr")){
0405 if(hist)hist->GetXaxis()->SetNdivisions(506);
0406 if(histZeroApe)histZeroApe->GetXaxis()->SetNdivisions(506);
0407 if(designHist)designHist->GetXaxis()->SetNdivisions(506);
0408 }
0409 if(sectorName.Contains("h_hitsPixel")){
0410 TGaxis::SetMaxDigits(3);
0411 }
0412
0413 this->setLineWidth(v_hist, 2);
0414 if(histZeroApe)histZeroApe->SetLineColor(2);
0415 if(histZeroApe)histZeroApe->SetLineStyle(2);
0416 if(hist)hist->SetLineColor(4);
0417 if(hist)hist->SetLineStyle(4);
0418
0419 TCanvas* canvas = new TCanvas("canvas");
0420 canvas->cd();
0421
0422 this->draw(v_hist);
0423 if(delta0_){
0424 const double xLow(designHist->GetXaxis()->GetXmin());
0425 const double xUp(designHist->GetXaxis()->GetXmax());
0426 TLine* baseline(nullptr);
0427 baseline = new TLine(xLow,*delta0_,xUp,*delta0_);
0428 baseline->SetLineStyle(2);
0429 baseline->SetLineWidth(2);
0430
0431 baseline->DrawLine(xLow,*delta0_,xUp,*delta0_);
0432 delete baseline;
0433 }
0434 canvas->Modified();
0435 canvas->Update();
0436
0437
0438 TPaveStats* statsZeroApe(nullptr);
0439 if(histZeroApe)statsZeroApe = (TPaveStats*)histZeroApe->GetListOfFunctions()->FindObject("stats");
0440 if(statsZeroApe){
0441 statsZeroApe->SetY1NDC(0.58);
0442 statsZeroApe->SetY2NDC(0.78);
0443 statsZeroApe->SetLineColor(2);
0444 statsZeroApe->SetTextColor(2);
0445 }
0446 TPaveStats* stats(nullptr);
0447 if(hist)stats = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
0448 if(stats){
0449 stats->SetY1NDC(0.37);
0450 stats->SetY2NDC(0.57);
0451 stats->SetLineColor(4);
0452 stats->SetTextColor(4);
0453 }
0454 canvas->Modified();
0455 canvas->Update();
0456
0457 TLegend* legend(nullptr);
0458 const LegendEntries& legendEntries = this->adjustLegendEntry(sectorName, hist, histZeroApe, designHist);
0459 legend = new TLegend(legendXmin_, legendYmin_, legendXmax_, legendYmax_);
0460 this->adjustLegend(legend);
0461 legend->AddEntry(designHist, legendEntries.designLegendEntry, "l");
0462 if(histZeroApe)legend->AddEntry(histZeroApe, legendEntries.legendEntryZeroApe, "l");
0463 if(hist)legend->AddEntry(hist, legendEntries.legendEntry, "l");
0464 legend->Draw("same");
0465
0466 canvas->Modified();
0467 canvas->Update();
0468
0469 const TString plotName(outpath_->Copy().Append(sectorName));
0470 canvas->Print(plotName + ".eps");
0471 canvas->Print(plotName + ".png");
0472
0473 if(legend)legend->Delete();
0474 if(stats)stats->Delete();
0475 if(canvas)canvas->Close();
0476 this->cleanup(v_hist);
0477 if(designHist)designHist->Delete();
0478 if(histZeroApe)histZeroApe->Delete();
0479 if(hist)hist->Delete();
0480 }
0481
0482
0483 void DrawPlot::adjustLegend(TLegend*& legend)const{
0484 if(!thesisMode_){
0485 legend->SetFillColor(0);
0486 legend->SetFillStyle(0);
0487 legend->SetTextSize(0.04);
0488 legend->SetMargin(0.12);
0489 legend->SetFillStyle(1001);
0490
0491 }
0492 else{
0493 legend->SetFillColor(0);
0494 legend->SetFillStyle(0);
0495 legend->SetTextSize(0.035);
0496 legend->SetMargin(0.12);
0497 legend->SetFillStyle(1001);
0498 }
0499 }
0500
0501
0502
0503
0504
0505 void DrawPlot::scale(std::vector<TH1*>& v_hist, const double factor)const{
0506 for(std::vector<TH1*>::const_iterator i_hist = v_hist.begin(); i_hist != v_hist.end(); ++i_hist){
0507 TH1* hist(*i_hist);
0508 const double integral(hist->Integral(0,hist->GetNbinsX()+1));
0509 if(integral>0.)hist->Scale(factor/integral);
0510 hist->SetYTitle(TString(hist->GetYaxis()->GetTitle()) + " (scaled)");
0511 }
0512 }
0513
0514
0515 double DrawPlot::maximumY(std::vector<TH1*>& v_hist)const{
0516 double maxY(-999.);
0517 for(std::vector<TH1*>::const_iterator i_hist = v_hist.begin(); i_hist != v_hist.end(); ++i_hist){
0518 const TH1* hist(*i_hist);
0519 const double histY(hist->GetMaximum());
0520 if(i_hist==v_hist.begin())maxY = histY;
0521 else if(maxY<histY)maxY = histY;
0522 }
0523 return maxY;
0524 }
0525
0526
0527 double DrawPlot::minimumY(std::vector<TH1*>& v_hist)const{
0528 double minY(-999.);
0529 for(std::vector<TH1*>::const_iterator i_hist = v_hist.begin(); i_hist != v_hist.end(); ++i_hist){
0530 const TH1* hist(*i_hist);
0531 const double histY(hist->GetMinimum());
0532 if(i_hist==v_hist.begin())minY = histY;
0533 else if(minY<histY)minY = histY;
0534 }
0535 return minY;
0536 }
0537
0538
0539 void DrawPlot::setRangeUser(std::vector<TH1*>& v_hist, const double minY, const double maxY)const{
0540 for(std::vector<TH1*>::const_iterator i_hist = v_hist.begin(); i_hist != v_hist.end(); ++i_hist){
0541 TH1* hist(*i_hist);
0542 hist->GetYaxis()->SetRangeUser(minY, maxY);
0543 }
0544 }
0545
0546
0547 void DrawPlot::setLineWidth(std::vector<TH1*>& v_hist, const unsigned int width)const{
0548 for(std::vector<TH1*>::const_iterator i_hist = v_hist.begin(); i_hist != v_hist.end(); ++i_hist){
0549 TH1* hist(*i_hist);
0550 hist->SetLineWidth(width);
0551 }
0552 }
0553
0554
0555 void DrawPlot::draw(std::vector<TH1*>& v_hist)const{
0556 for(std::vector<TH1*>::const_iterator i_hist = v_hist.begin(); i_hist != v_hist.end(); ++i_hist){
0557 TH1* hist(*i_hist);
0558 if(i_hist==v_hist.begin())hist->Draw();
0559 else hist->Draw("sameS");
0560 }
0561 }
0562
0563
0564 void DrawPlot::cleanup(std::vector<TH1*>& v_hist)const{
0565 for(std::vector<TH1*>::iterator i_hist = v_hist.begin(); i_hist != v_hist.end(); ++i_hist){
0566 TH1* hist(*i_hist);
0567 hist = 0;
0568 }
0569 v_hist.clear();
0570 }
0571
0572
0573
