Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:44:20

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   //if(!file_ || !fileZeroApe_ || !designFile_ || !baselineFile){
0055     // Not needed: root gives error by default when file is not found
0056     //std::cout<<"\n\tInput file not found, please check file name: "<<*fileName<<"\n";
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     //mode = "";
0181     //precision = ;
0182     //unit = "";
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     //std::cout<<"Sector: "<<ss_sector.str()<<"\n";
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   //const double minY(this->minimumY(v_hist));
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     //baseline = new TLine(0.0005,*delta0_,0.01,*delta0_); baseline->Draw("same");
0431     baseline->DrawLine(xLow,*delta0_,xUp,*delta0_);
0432     delete baseline;
0433   }
0434   canvas->Modified();
0435   canvas->Update();
0436   
0437   //TPaveStats* stats =(TPaveStats*) hist->GetListOfFunctions()->At(1);
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     //legend->SetBorderSize(0);
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