Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-28 03:58:37

0001 #include "Alignment/OfflineValidation/macros/PlotAlignmentValidation.h"
0002 
0003 #include "Alignment/OfflineValidation/macros/TkAlStyle.cc"
0004 #include "Alignment/OfflineValidation/interface/TkOffTreeVariables.h"
0005 
0006 #include "Math/ProbFunc.h"
0007 
0008 #include "TAxis.h"
0009 #include "TCanvas.h"
0010 #include "TDirectory.h"
0011 #include "TDirectoryFile.h"
0012 #include "TF1.h"
0013 #include "TFile.h"
0014 #include "TGaxis.h"
0015 #include "TH2F.h"
0016 #include "THStack.h"
0017 #include "TKey.h"
0018 #include "TLatex.h"
0019 #include "TLegend.h"
0020 #include "TLegendEntry.h"
0021 #include "TPad.h"
0022 #include "TPaveStats.h"
0023 #include "TPaveText.h"
0024 #include "TProfile.h"
0025 #include "TRandom3.h"
0026 #include "TRegexp.h"
0027 #include "TROOT.h"
0028 #include "TString.h"
0029 #include "TStyle.h"
0030 #include "TSystem.h"
0031 #include "TTree.h"
0032 
0033 #include <algorithm>
0034 #include <cmath>
0035 #include <cstdio>
0036 #include <cstdlib>
0037 #include <exception>
0038 #include <fstream>
0039 #include <iostream>
0040 #include <memory>
0041 #include <sstream>
0042 #include <string>
0043 #include <vector>
0044 
0045 /*! \class PlotAlignmentValidation
0046  *  \brief Class PlotAlignmentValidation
0047  *         Class used as the last step for Offline Track Validation tool.
0048  *         The main goal of this class is creating the plots regarding DMRs and Surface Deformations for modules and substructures.
0049  */
0050 
0051 //------------------------------------------------------------------------------
0052 /*! \fn PlotAlignmentValidation
0053  *  \brief Constructor for the class
0054  */
0055 
0056 PlotAlignmentValidation::PlotAlignmentValidation(bool bigtext) : bigtext_(bigtext) {
0057   setOutputDir(".");
0058   setTreeBaseDir();
0059   sourcelist = NULL;
0060 
0061   moreThanOneSource = false;
0062   useFit_ = false;
0063 
0064   // Force ROOT to use scientific notation even with smaller datasets
0065   TGaxis::SetMaxDigits(4);
0066   // (This sets a static variable: correct in .eps images but must be set
0067   // again manually when viewing the .root files)
0068 
0069   // Make ROOT calculate histogram statistics using all data,
0070   // regardless of displayed range
0071   TH1::StatOverflows(kTRUE);
0072 
0073   //show all information in the legend by default
0074   legendOptions(TkAlStyle::legendoptions);
0075 }
0076 
0077 //------------------------------------------------------------------------------
0078 /*! \fn PlotAlignmentValidation
0079  *  \brief Constructor for the class. This function also retrieves the list of root files used to produce DMRs and Surface Deformations
0080  */
0081 PlotAlignmentValidation::PlotAlignmentValidation(
0082     const char* inputFile, std::string legendName, int lineColor, int lineStyle, bool bigtext)
0083     : PlotAlignmentValidation(bigtext) {
0084   loadFileList(inputFile, legendName, lineColor, lineStyle);
0085 }
0086 
0087 //------------------------------------------------------------------------------
0088 /*! \fn ~PlotAlignmentValidation
0089  *  \brief Default destructor
0090  */
0091 PlotAlignmentValidation::~PlotAlignmentValidation() {
0092   for (std::vector<TkOfflineVariables*>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
0093     delete (*it);
0094   }
0095 
0096   delete sourcelist;
0097 }
0098 
0099 //------------------------------------------------------------------------------
0100 /*!
0101  * \fn openSummaryFile
0102  * \brief Create/open the root and txt summary files, where the DMR histograms and the associtated mean/sigma are stored respectively
0103  */
0104 
0105 void PlotAlignmentValidation::openSummaryFile() {
0106   if (!openedsummaryfile) {
0107     openedsummaryfile = true;
0108     summaryfile.open(TString(outputDir + "/") + summaryfilename + ".txt");
0109     //Rootfile introduced to store the DMR histograms
0110     rootsummaryfile = new TFile(TString(outputDir + "/") + summaryfilename + ".root", "RECREATE");
0111 
0112     for (auto vars : sourceList) {
0113       summaryfile << "\t" << vars->getName();
0114     }
0115     summaryfile << "\tformat={}\n";
0116   } else {
0117     //Check for the rootfile to be open, and open it in case it is not already.
0118     if (!rootsummaryfile->IsOpen())
0119       rootsummaryfile->Open(TString(outputDir + "/") + summaryfilename + ".root", "UPDATE");
0120   }
0121 }
0122 
0123 //------------------------------------------------------------------------------
0124 /*! \fn loadFileList
0125  *  \brief Add to the list of sources the rootfile associated to a particular geometry
0126  */
0127 void PlotAlignmentValidation::loadFileList(const char* inputFile, std::string legendName, int lineColor, int lineStyle) {
0128   if (openedsummaryfile) {
0129     std::cout << "Can't load a root file after opening the summary file!" << std::endl;
0130     assert(0);
0131   }
0132   sourceList.push_back(new TkOfflineVariables(inputFile, treeBaseDir, legendName, lineColor, lineStyle));
0133 }
0134 
0135 //------------------------------------------------------------------------------
0136 /*! \fn useFitForDMRplots
0137  *  \brief Store the selected boolean in one of the private members of the class
0138  */
0139 void PlotAlignmentValidation::useFitForDMRplots(bool usefit) { useFit_ = usefit; }
0140 
0141 //------------------------------------------------------------------------------
0142 /*! \fn numberOfLayers
0143  *  \brief Select the number of layers associated to a subdetector.
0144  */
0145 //TODO Possible improvement: reduce the number of switches in the code by implementing a map
0146 int PlotAlignmentValidation::numberOfLayers(int phase, int subdetector) {
0147   switch (phase) {
0148     case 0:
0149       switch (subdetector) {
0150         case 1:
0151           return 3;
0152         case 2:
0153           return 2;
0154         case 3:
0155           return 4;
0156         case 4:
0157           return 3;
0158         case 5:
0159           return 6;
0160         case 6:
0161           return 9;
0162         default:
0163           assert(false);
0164       }
0165     case 1:
0166       switch (subdetector) {
0167         case 1:
0168           return 4;
0169         case 2:
0170           return 3;
0171         case 3:
0172           return 4;
0173         case 4:
0174           return 3;
0175         case 5:
0176           return 6;
0177         case 6:
0178           return 9;
0179         default:
0180           assert(false);
0181       }
0182     default:
0183       assert(false);
0184   }
0185   return 0;
0186 }
0187 
0188 //------------------------------------------------------------------------------
0189 /*! \fn maxNumberOfLayers
0190  *  \brief Return the number of layers of a subdetector
0191  */
0192 int PlotAlignmentValidation::maxNumberOfLayers(int subdetector) {
0193   int result = 0;
0194   for (auto it = sourceList.begin(); it != sourceList.end(); ++it) {
0195     result = max(result, numberOfLayers((*it)->getPhase(), subdetector));
0196   }
0197   return result;
0198 }
0199 
0200 //------------------------------------------------------------------------------
0201 /*! \fn legendOptions
0202  *  \brief Assign legend options to members of the class
0203  */
0204 void PlotAlignmentValidation::legendOptions(TString options) {
0205   showMean_ = false;
0206   showRMS_ = false;
0207   showMeanError_ = false;
0208   showRMSError_ = false;
0209   showModules_ = false;
0210   showUnderOverFlow_ = false;
0211   options.ReplaceAll(" ", "").ToLower();
0212   if (options.Contains("mean") || options.Contains("all"))
0213     showMean_ = true;
0214   if (options.Contains("meanerror") || options.Contains("all"))
0215     showMeanError_ = true;
0216   if (options.Contains("rms") || options.Contains("all"))
0217     showRMS_ = true;
0218   if (options.Contains("rmserror") || options.Contains("all"))
0219     showRMSError_ = true;
0220   if (options.Contains("modules") || options.Contains("all"))
0221     showModules_ = true;
0222   if (options.Contains("under") || options.Contains("over") || options.Contains("outside") || options.Contains("all"))
0223     showUnderOverFlow_ = true;
0224 
0225   twolines_ = (showUnderOverFlow_ && (showMean_ + showMeanError_ + showRMS_ + showRMSError_ >= 1) && bigtext_);
0226 }
0227 
0228 //------------------------------------------------------------------------------
0229 /*! \fn setOutputDir
0230  *  \brief Set the output direcotry
0231  */
0232 void PlotAlignmentValidation::setOutputDir(std::string dir) {
0233   if (openedsummaryfile) {
0234     std::cout << "Can't set the output dir after opening the summary file!" << std::endl;
0235     assert(0);
0236   }
0237   outputDir = dir;
0238   gSystem->mkdir(outputDir.data(), true);
0239 }
0240 
0241 //------------------------------------------------------------------------------
0242 /*! \fn plotSubDetResiduals
0243  *  \brief Function used to plot residuals for a subdetector
0244  */
0245 void PlotAlignmentValidation::plotSubDetResiduals(bool plotNormHisto, unsigned int subDetId) {
0246   gStyle->SetOptStat(11111);
0247   gStyle->SetOptFit(0000);
0248 
0249   TCanvas* c = new TCanvas("c", "c");
0250   c->SetTopMargin(0.15);
0251   TString histoName = "";
0252   if (plotNormHisto) {
0253     histoName = "h_NormXprime";
0254   } else
0255     histoName = "h_Xprime_";
0256   switch (subDetId) {
0257     case 1:
0258       histoName += "TPBBarrel_0";
0259       break;
0260     case 2:
0261       histoName += "TPEendcap_1";
0262       break;
0263     case 3:
0264       histoName += "TPEendcap_2";
0265       break;
0266     case 4:
0267       histoName += "TIBBarrel_0";
0268       break;
0269     case 5:
0270       histoName += "TIDEndcap_1";
0271       break;
0272     case 6:
0273       histoName += "TIDEndcap_2";
0274       break;
0275     case 7:
0276       histoName += "TOBBarrel_3";
0277       break;
0278     case 8:
0279       histoName += "TECEndcap_4";
0280       break;
0281     case 9:
0282       histoName += "TECEndcap_5";
0283       break;
0284   }
0285   int tmpcounter = 0;
0286   TH1* sumHisto = 0;
0287   for (std::vector<TkOfflineVariables*>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
0288     if (tmpcounter == 0) {
0289       TFile* f = (*it)->getFile();
0290       sumHisto = (TH1*)f->FindKeyAny(histoName)->ReadObj();  //FindObjectAny(histoName.Data());
0291       sumHisto->SetLineColor(tmpcounter + 1);
0292       sumHisto->SetLineStyle(tmpcounter + 1);
0293       sumHisto->GetFunction("tmp")->SetBit(TF1::kNotDraw);
0294       sumHisto->Draw();
0295 
0296       //get statistic box coordinate to plot all boxes one below the other
0297       //gStyle->SetStatY(0.91);
0298       //gStyle->SetStatW(0.15);
0299       //gStyle->SetStatBorderSize(1);
0300       //gStyle->SetStatH(0.10);
0301 
0302       tmpcounter++;
0303     } else {
0304       sumHisto = (TH1*)(*it)->getFile()->FindObjectAny(histoName);
0305       sumHisto->SetLineColor(tmpcounter + 1);
0306       sumHisto->SetLineStyle(tmpcounter + 1);
0307       sumHisto->GetFunction("tmp")->SetBit(TF1::kNotDraw);
0308       //hstack->Add(sumHisto);
0309 
0310       c->Update();
0311       tmpcounter++;
0312     }
0313     TObject* statObj = sumHisto->GetListOfFunctions()->FindObject("stats");
0314     if (statObj && statObj->InheritsFrom(TPaveStats::Class())) {
0315       TPaveStats* stats = static_cast<TPaveStats*>(statObj);
0316       stats->SetLineColor(tmpcounter + 1);
0317       stats->SetTextColor(tmpcounter + 1);
0318       stats->SetFillColor(10);
0319       stats->SetX1NDC(0.91 - tmpcounter * 0.1);
0320       stats->SetX2NDC(0.15);
0321       stats->SetY1NDC(1);
0322       stats->SetY2NDC(0.10);
0323       sumHisto->Draw("sames");
0324     }
0325   }
0326   //hstack->Draw("nostack");
0327   char PlotName[1000];
0328   sprintf(PlotName, "%s/%s.png", outputDir.c_str(), histoName.Data());
0329   c->Print(PlotName);
0330   sprintf(PlotName, "%s/%s.eps", outputDir.c_str(), histoName.Data());
0331   c->Print(PlotName);
0332   sprintf(PlotName, "%s/%s.pdf", outputDir.c_str(), histoName.Data());
0333   c->Print(PlotName);
0334   sprintf(PlotName, "%s/%s.root", outputDir.c_str(), histoName.Data());
0335   c->Print(PlotName);
0336   //delete c;
0337   //c=0;
0338 }
0339 
0340 //------------------------------------------------------------------------------
0341 void PlotAlignmentValidation::plotHitMaps() {
0342   //gStyle->SetOptStat(0);
0343 
0344   TCanvas* c = new TCanvas("c", "c", 1200, 400);
0345   c->Divide(3, 1);
0346   //ps->NewPage();
0347 
0348   //-------------------------------------------------
0349   //plot Hit map
0350   //-------------------------------------------------
0351   std::string histName_ = "Entriesprofile";
0352   c->cd(1);
0353   TTree* tree = (*sourceList.begin())->getTree();
0354   tree->Draw("entries:posR:posZ", "", "COLZ2Prof");
0355   c->cd(2);
0356   tree->Draw("entries:posY:posX", "", "COLZ2Prof");
0357   c->cd(3);
0358   tree->Draw("entries:posR:posPhi", "", "COLZ2Prof");
0359 
0360   char PlotName[1000];
0361   sprintf(PlotName, "%s/%s.png", outputDir.c_str(), histName_.c_str());
0362   c->Print(PlotName);
0363   sprintf(PlotName, "%s/%s.eps", outputDir.c_str(), histName_.c_str());
0364   c->Print(PlotName);
0365   sprintf(PlotName, "%s/%s.pdf", outputDir.c_str(), histName_.c_str());
0366   c->Print(PlotName);
0367   sprintf(PlotName, "%s/%s.root", outputDir.c_str(), histName_.c_str());
0368   c->Print(PlotName);
0369   //   //c->Update();
0370   c->Close();
0371   //----------------------------------------------------
0372 }
0373 
0374 //------------------------------------------------------------------------------
0375 void PlotAlignmentValidation::plotOutlierModules(const char* outputFileName,
0376                                                  std::string plotVariable,
0377                                                  float plotVariable_cut,
0378                                                  int unsigned minHits) {
0379   Int_t counter = 0;
0380 
0381   gStyle->SetOptStat(111111);
0382   gStyle->SetStatY(0.9);
0383   //TList* treelist=getTreeList();
0384 
0385   TCanvas* c1 = new TCanvas("canv", "canv", 800, 500);
0386   outputFile = outputDir + '/' + outputFileName;
0387   c1->Print((outputFile + '[').Data());
0388 
0389   c1->Divide(2, 1);
0390 
0391   TTree* tree = (*sourceList.begin())->getTree();
0392   TkOffTreeVariables* treeMem = 0;  // ROOT will initilise
0393   tree->SetBranchAddress("TkOffTreeVariables", &treeMem);
0394 
0395   Long64_t nentries = tree->GetEntriesFast();
0396 
0397   for (Long64_t i = 0; i < nentries; i++) {
0398     tree->GetEntry(i);
0399     float var = 0;
0400     if (plotVariable == "chi2PerDofX")
0401       var = treeMem->chi2PerDofX;
0402     else if (plotVariable == "chi2PerDofY")
0403       var = treeMem->chi2PerDofY;
0404     else if (plotVariable == "fitMeanX")
0405       var = treeMem->fitMeanX;
0406     else if (plotVariable == "fitMeanY")
0407       var = treeMem->fitMeanY;
0408     else if (plotVariable == "fitSigmaX")
0409       var = treeMem->fitSigmaX;
0410     else if (plotVariable == "fitSigmaY")
0411       var = treeMem->fitSigmaY;
0412     else {
0413       cout << "There is no variable " << plotVariable << " included in the tree." << endl;
0414       break;
0415     }
0416     //   cout<<"treeMem->entries  "<<treeMem->entries<<endl;
0417     //  cout<<"var                  "<<var<<endl;
0418     //  cout<<"plotVariable_cut     "<<plotVariable_cut<<endl;
0419 
0420     if (var > plotVariable_cut && treeMem->entries > minHits) {
0421       TFile* f = (*sourceList.begin())->getFile();  //(TFile*)sourcelist->First();
0422 
0423       if (f->FindKeyAny(treeMem->histNameX.c_str()) != 0) {
0424         TH1* h =
0425             (TH1*)f->FindKeyAny(treeMem->histNameX.c_str())->ReadObj();  //f->FindObjectAny(treeMem->histNameX.c_str());
0426         gStyle->SetOptFit(0111);
0427         cout << "hist name " << h->GetName() << endl;
0428 
0429         TString path = (char*)strstr(gDirectory->GetPath(), "TrackerOfflineValidation");
0430         //cout<<"hist path "<<path<<endl;
0431         //cout<<"wrote text "<<endl;
0432         if (h)
0433           cout << h->GetEntries() << endl;
0434 
0435         //modules' location as title
0436         c1->cd(0);
0437         TPaveText* text = new TPaveText(0, 0.95, 0.99, 0.99);
0438         text->AddText(path);
0439         text->SetFillColor(0);
0440         text->SetShadowColor(0);
0441         text->SetBorderSize(0);
0442         text->Draw();
0443 
0444         //residual histogram
0445         c1->cd(1);
0446         TPad* subpad = (TPad*)c1->GetPad(1);
0447         subpad->SetPad(0, 0, 0.5, 0.94);
0448         h->Draw();
0449 
0450         //norm. residual histogram
0451         h = (TH1*)f->FindObjectAny(treeMem->histNameNormX.c_str());
0452         if (h)
0453           cout << h->GetEntries() << endl;
0454         c1->cd(2);
0455         TPad* subpad2 = (TPad*)c1->GetPad(2);
0456         subpad2->SetPad(0.5, 0, 0.99, 0.94);
0457         h->Draw();
0458 
0459         c1->Print(outputFile);
0460         counter++;
0461       } else {
0462         cout << "There are no residual histograms on module level stored!" << endl;
0463         cout << "Please make sure that moduleLevelHistsTransient = cms.bool(False) in the validation job!" << endl;
0464         break;
0465       }
0466     }
0467   }
0468   c1->Print((outputFile + "]").Data());
0469   if (counter == 0)
0470     cout << "no bad modules found" << endl;
0471 
0472   //read the number of entries in the t3
0473   //TTree* tree=0;
0474   //tree=(TTree*)treeList->At(0);
0475 
0476   //c1->Close();
0477 }
0478 
0479 //------------------------------------------------------------------------------
0480 /*! \fn getTreeList
0481  *  \brief Extract from the rootfiles stored in the sourcelist the corresponding trees.
0482  */
0483 TList* PlotAlignmentValidation::getTreeList() {
0484   TList* treeList = new TList();
0485   TFile* first_source = (TFile*)sourcelist->First();
0486   std::cout << first_source->GetName() << std::endl;
0487   TDirectoryFile* d = (TDirectoryFile*)first_source->Get(treeBaseDir.c_str());
0488   treeList->Add((TTree*)(*d).Get("TkOffVal"));
0489 
0490   if (moreThanOneSource == true) {
0491     TFile* nextsource = (TFile*)sourcelist->After(first_source);
0492     while (nextsource) {
0493       std::cout << nextsource->GetName() << std::endl;
0494       d = (TDirectoryFile*)nextsource->Get("TrackerOfflineValidation");
0495 
0496       treeList->Add((TTree*)(*d).Get("TkOffVal"));
0497 
0498       nextsource = (TFile*)sourcelist->After(nextsource);
0499     }
0500   }
0501   return treeList;
0502 }
0503 
0504 //------------------------------------------------------------------------------
0505 void PlotAlignmentValidation::setTreeBaseDir(std::string dir) { treeBaseDir = dir; }
0506 
0507 //------------------------------------------------------------------------------
0508 void PlotAlignmentValidation::plotSurfaceShapes(const std::string& options, const std::string& residType) {
0509   cout << "-------- plotSurfaceShapes called with " << options << endl;
0510   if (options == "none")
0511     return;
0512   else if (options == "coarse") {
0513     plotSS("subdet=1");
0514     plotSS("subdet=2");
0515     plotSS("subdet=3");
0516     plotSS("subdet=4");
0517     plotSS("subdet=5");
0518     plotSS("subdet=6");
0519   }
0520   // else if (options == "fine") ...
0521   else
0522     plotSS(options, residType);
0523 
0524   return;
0525 }
0526 
0527 //------------------------------------------------------------------------------
0528 void PlotAlignmentValidation::plotSS(const std::string& options, const std::string& residType) {
0529   if (residType == "") {
0530     plotSS(options, "ResXvsXProfile");
0531     plotSS(options, "ResXvsYProfile");
0532     return;
0533   }
0534 
0535   int bkperrorx = gStyle->GetErrorX();
0536   gStyle->SetErrorX(1);  //regardless of style settings, we want x error bars here
0537 
0538   int plotLayerN = 0;
0539   //  int plotRingN  = 0;
0540   //  bool plotPlain = false;
0541   bool plotLayers = false;  // overrides plotLayerN
0542   //  bool plotRings  = false;  // Todo: implement this?
0543   bool plotSplits = false;
0544   int plotSubDetN = 0;  // if zero, plot all
0545 
0546   TRegexp layer_re("layer=[0-9]+");
0547   Ssiz_t index, len;
0548   if (options.find("split") != std::string::npos) {
0549     plotSplits = true;
0550   }
0551   if (options.find("layers") != std::string::npos) {
0552     plotLayers = true;
0553   }
0554   if ((index = layer_re.Index(options, &len)) != -1) {
0555     if (plotLayers) {
0556       std::cerr << "Warning: option 'layers' overrides 'layer=N'" << std::endl;
0557     } else {
0558       std::string substr = options.substr(index + 6, len - 6);
0559       plotLayerN = atoi(substr.c_str());
0560     }
0561   }
0562 
0563   TRegexp subdet_re("subdet=[1-6]+");
0564   if ((index = subdet_re.Index(options, &len)) != -1) {
0565     std::string substr = options.substr(index + 7, len - 7);
0566     plotSubDetN = atoi(substr.c_str());
0567   }
0568 
0569   gStyle->SetOptStat(0);
0570 
0571   TCanvas c("canv", "canv");
0572 
0573   // todo: title, min/max, nbins?
0574 
0575   // Loop over detectors
0576   for (int iSubDet = 1; iSubDet <= 6; ++iSubDet) {
0577     // TEC requires special care since rings 1-4 and 5-7 are plotted separately
0578     bool isTEC = (iSubDet == 6);
0579 
0580     // if subdet is specified, skip other subdets
0581     if (plotSubDetN != 0 && iSubDet != plotSubDetN)
0582       continue;
0583 
0584     // Skips plotting too high layers
0585     // if it's a mixture of phase 0 and 1, the phase 0 files will be skipped
0586     //  when plotting the higher layers of BPIX and FPIX
0587     if (plotLayerN > maxNumberOfLayers(iSubDet)) {
0588       continue;
0589     }
0590 
0591     int minlayer = plotLayers ? 1 : plotLayerN;
0592     int maxlayer = plotLayers ? maxNumberOfLayers(iSubDet) : plotLayerN;
0593     // see later where this is used
0594     int maxlayerphase0 = plotLayers ? numberOfLayers(0, iSubDet) : plotLayerN;
0595 
0596     for (int layer = minlayer; layer <= maxlayer; layer++) {
0597       // two plots for TEC, skip first
0598       for (int iTEC = 0; iTEC < 2; iTEC++) {
0599         if (!isTEC && iTEC == 0)
0600           continue;
0601 
0602         char selection[1000];
0603         if (!isTEC) {
0604           if (layer == 0)
0605             sprintf(selection, "subDetId==%d", iSubDet);
0606           else
0607             sprintf(selection, "subDetId==%d && layer == %d", iSubDet, layer);
0608         } else {          // TEC
0609           if (iTEC == 0)  // rings
0610             sprintf(selection, "subDetId==%d && ring <= 4", iSubDet);
0611           else
0612             sprintf(selection, "subDetId==%d && ring > 4", iSubDet);
0613         }
0614 
0615         // Title for plot and name for the file
0616 
0617         TString subDetName;
0618         switch (iSubDet) {
0619           case 1:
0620             subDetName = "BPIX";
0621             break;
0622           case 2:
0623             subDetName = "FPIX";
0624             break;
0625           case 3:
0626             subDetName = "TIB";
0627             break;
0628           case 4:
0629             subDetName = "TID";
0630             break;
0631           case 5:
0632             subDetName = "TOB";
0633             break;
0634           case 6:
0635             subDetName = "TEC";
0636             break;
0637         }
0638 
0639         TString secondline = "";
0640         if (layer != 0) {
0641           // TEC and TID have discs, the rest have layers
0642           if (iSubDet == 4 || iSubDet == 6)
0643             secondline = "disc ";
0644           else {
0645             secondline = "layer ";
0646           }
0647           secondline += Form("%d", layer);
0648           secondline += " ";
0649         }
0650         if (isTEC && iTEC == 0)
0651           secondline += TString("R1-4");
0652         if (isTEC && iTEC > 0)
0653           secondline += TString("R5-7");
0654 
0655         // Generate histograms with selection
0656         TLegend* legend = 0;
0657         // Any file from phase 0 will be skipped if the last argument is false
0658         THStack* hs = addHists(selection, residType, &legend, false, /*validforphase0 = */ layer <= maxlayerphase0);
0659         if (!hs || hs->GetHists() == 0 || hs->GetHists()->GetSize() == 0) {
0660           std::cout << "No histogram for " << subDetName << ", perhaps not enough data? Creating default histogram."
0661                     << std::endl;
0662           if (hs == 0)
0663             hs = new THStack("hstack", "");
0664 
0665           TProfile* defhist = new TProfile("defhist", "Empty default histogram", 100, -1, 1, -1, 1);
0666           hs->Add(defhist);
0667           hs->Draw();
0668         } else {
0669           hs->Draw("nostack PE");
0670           modifySSHistAndLegend(hs, legend);
0671           legend->Draw();
0672           setTitleStyle(*hs, "", "", iSubDet, true, secondline);
0673 
0674           // Adjust Labels
0675           TH1* firstHisto = (TH1*)hs->GetHists()->First();
0676           TString xName = firstHisto->GetXaxis()->GetTitle();
0677           TString yName = firstHisto->GetYaxis()->GetTitle();
0678           hs->GetHistogram()->GetXaxis()->SetTitleColor(kBlack);
0679           hs->GetHistogram()->GetXaxis()->SetTitle(xName);
0680           hs->GetHistogram()->GetYaxis()->SetTitleColor(kBlack);
0681           // micrometers:
0682           yName.ReplaceAll("cm", "#mum");
0683           hs->GetHistogram()->GetYaxis()->SetTitle(yName);
0684         }
0685 
0686         // Save plot to file
0687         std::ostringstream plotName;
0688         plotName << outputDir << "/SurfaceShape_" << subDetName << "_";
0689         plotName << residType;
0690         if (layer != 0) {
0691           plotName << "_";
0692           // TEC and TID have discs, the rest have layers
0693           if (iSubDet == 4 || iSubDet == 6)
0694             plotName << "disc";
0695           else {
0696             plotName << "layer";
0697           }
0698           plotName << layer;
0699         }
0700         if (isTEC && iTEC == 0)
0701           plotName << "_"
0702                    << "R1-4";
0703         if (isTEC && iTEC > 0)
0704           plotName << "_"
0705                    << "R5-7";
0706 
0707         // PNG,EPS,PDF files
0708         c.Update();
0709         c.Print((plotName.str() + ".png").c_str());
0710         c.Print((plotName.str() + ".eps").c_str());
0711         c.Print((plotName.str() + ".pdf").c_str());
0712 
0713         // ROOT file
0714         TFile f((plotName.str() + ".root").c_str(), "recreate");
0715         c.Write();
0716         f.Close();
0717 
0718         delete legend;
0719         delete hs;
0720       }
0721     }
0722   }
0723   gStyle->SetErrorX(bkperrorx);
0724 
0725   return;
0726 }
0727 
0728 //------------------------------------------------------------------------------
0729 /*! \fn plotDMR
0730  *  \brief Main function used to plot DMRs for a single IOV printing the canvases in the output directory and saving histograms and fit funtions in a root file.
0731  */
0732 void PlotAlignmentValidation::plotDMR(const std::string& variable, Int_t minHits, const std::string& options) {
0733   // If several, comma-separated values are given in 'variable',
0734   // call plotDMR with each value separately.
0735   // If a comma is found, the string is divided to two.
0736   // (no space allowed)
0737   std::size_t findres = variable.find(",");
0738   if (findres != std::string::npos) {
0739     std::string substring1 = variable.substr(0, findres);
0740     std::string substring2 = variable.substr(findres + 1, std::string::npos);
0741     plotDMR(substring1, minHits, options);
0742     plotDMR(substring2, minHits, options);
0743     return;
0744   }
0745 
0746   // Variable name should end with X or Y. If it doesn't, recursively calls plotDMR twice with
0747   // X and Y added, respectively
0748   if (variable == "mean" || variable == "median" || variable == "meanNorm" || variable == "rms" ||
0749       variable == "rmsNorm") {
0750     plotDMR(variable + "X", minHits, options);
0751     plotDMR(variable + "Y", minHits, options);
0752     return;
0753   }
0754 
0755   // options:
0756   // -plain (default, the whole distribution)
0757   // -split (distribution splitted to two)
0758   // -layers (plain db for each layer/disc superimposed in one plot)
0759   // -layersSeparate (plain db for each layer/disc in separate plots)
0760   // -layersSplit (splitted db for each layers/disc in one plot)
0761   // -layersSplitSeparate (splitted db, for each layers/disc in separate plots)
0762 
0763   TRegexp layer_re("layer=[0-9]+");
0764   bool plotPlain = false, plotSplits = false, plotLayers = false;
0765   int plotLayerN = 0;
0766   Ssiz_t index, len;
0767   if (options.find("plain") != std::string::npos) {
0768     plotPlain = true;
0769   }
0770   if (options.find("split") != std::string::npos) {
0771     plotSplits = true;
0772   }
0773   if (options.find("layers") != std::string::npos) {
0774     plotLayers = true;
0775   }
0776   if ((index = layer_re.Index(options, &len)) != -1) {
0777     if (plotLayers) {
0778       std::cerr << "Warning: option 'layers' overrides 'layer=N'" << std::endl;
0779     } else {
0780       std::string substr = options.substr(index + 6, len - 6);
0781       plotLayerN = atoi(substr.c_str());
0782     }
0783   }
0784 
0785   // Defaults to plotting only plain plot if empty (or invalid)
0786   // option string is given
0787   if (!plotPlain && !plotSplits) {
0788     plotPlain = true;
0789   }
0790 
0791   // This boolean array tells for which detector modules to plot split DMR plots
0792   // They are plotted for BPIX, FPIX, TIB and TOB
0793   static bool plotSplitsFor[6] = {true, true, true, false, true, false};
0794 
0795   DMRPlotInfo plotinfo;
0796 
0797   gStyle->SetOptStat(0);
0798 
0799   TCanvas c("canv", "canv");
0800 
0801   plotinfo.variable = variable;
0802   plotinfo.minHits = minHits;
0803   plotinfo.plotPlain = plotPlain;
0804   plotinfo.plotLayers = plotLayers;
0805 
0806   // width in cm
0807   // for DMRS, use 100 bins in range +-10 um, bin width 0.2um
0808   // if modified, check also TrackerOfflineValidationSummary_cfi.py and TrackerOfflineValidation_Standalone_cff.py
0809   if (variable == "meanX") {
0810     plotinfo.nbins = 50;
0811     plotinfo.min = -0.001;
0812     plotinfo.max = 0.001;
0813   } else if (variable == "meanY") {
0814     plotinfo.nbins = 50;
0815     plotinfo.min = -0.005;
0816     plotinfo.max = 0.005;
0817   } else if (variable == "medianX")
0818     if (plotSplits) {
0819       plotinfo.nbins = 50;
0820       plotinfo.min = -0.0005;
0821       plotinfo.max = 0.0005;
0822     } else {
0823       plotinfo.nbins = 100;
0824       plotinfo.min = -0.001;
0825       plotinfo.max = 0.001;
0826     }
0827   else if (variable == "medianY")
0828     if (plotSplits) {
0829       plotinfo.nbins = 50;
0830       plotinfo.min = -0.0005;
0831       plotinfo.max = 0.0005;
0832     } else {
0833       plotinfo.nbins = 100;
0834       plotinfo.min = -0.001;
0835       plotinfo.max = 0.001;
0836     }
0837   else if (variable == "meanNormX") {
0838     plotinfo.nbins = 100;
0839     plotinfo.min = -2.0;
0840     plotinfo.max = 2.0;
0841   } else if (variable == "meanNormY") {
0842     plotinfo.nbins = 100;
0843     plotinfo.min = -2.0;
0844     plotinfo.max = 2.0;
0845   } else if (variable == "rmsX") {
0846     plotinfo.nbins = 100;
0847     plotinfo.min = 0.0;
0848     plotinfo.max = 0.1;
0849   } else if (variable == "rmsY") {
0850     plotinfo.nbins = 100;
0851     plotinfo.min = 0.0;
0852     plotinfo.max = 0.1;
0853   } else if (variable == "rmsNormX") {
0854     plotinfo.nbins = 100;
0855     plotinfo.min = 0.3;
0856     plotinfo.max = 1.8;
0857   } else if (variable == "rmsNormY") {
0858     plotinfo.nbins = 100;
0859     plotinfo.min = 0.3;
0860     plotinfo.max = 1.8;
0861   } else {
0862     std::cerr << "Unknown variable " << variable << std::endl;
0863     plotinfo.nbins = 100;
0864     plotinfo.min = -0.1;
0865     plotinfo.max = 0.1;
0866   }
0867   //Begin loop on structures
0868   for (int i = 1; i <= 6; ++i) {
0869     // Skip strip detectors if plotting any "Y" variable
0870     if (i != 1 && i != 2 && variable.length() > 0 && variable[variable.length() - 1] == 'Y') {
0871       continue;
0872     }
0873 
0874     // Skips plotting too high layers
0875     if (plotLayerN > maxNumberOfLayers(i)) {
0876       continue;
0877     }
0878 
0879     plotinfo.plotSplits = plotSplits && plotSplitsFor[i - 1];
0880     if (!plotinfo.plotPlain && !plotinfo.plotSplits) {
0881       continue;
0882     }
0883 
0884     // Sets dimension of legend according to the number of plots
0885 
0886     bool hasheader = (TkAlStyle::legendheader != "");
0887 
0888     int nPlots = 1;
0889     if (plotinfo.plotSplits) {
0890       nPlots = 3;
0891     }
0892     // This will make the legend a bit bigger than necessary if there is a mixture of phase 0 and phase 1.
0893     // Not worth it to implement more complicated logic.
0894     if (plotinfo.plotLayers) {
0895       nPlots *= maxNumberOfLayers(i);
0896     }
0897     nPlots *= sourceList.size();
0898     if (twolines_) {
0899       nPlots *= 2;
0900     }
0901     nPlots += hasheader;
0902 
0903     double legendY = 0.80;
0904     if (nPlots > 3) {
0905       legendY -= 0.01 * (nPlots - 3);
0906     }
0907     if (bigtext_) {
0908       legendY -= 0.05;
0909     }
0910     if (legendY < 0.6) {
0911       std::cerr << "Warning: Huge legend!" << std::endl;
0912       legendY = 0.6;
0913     }
0914 
0915     THStack hstack("hstack", "hstack");
0916     plotinfo.maxY = 0;
0917     plotinfo.subDetId = i;
0918     plotinfo.legend = new TLegend(0.17, legendY, 0.85, 0.88);
0919     plotinfo.legend->SetNColumns(2);
0920     if (hasheader)
0921       plotinfo.legend->SetHeader(TkAlStyle::legendheader);
0922     if (bigtext_)
0923       plotinfo.legend->SetTextSize(TkAlStyle::textSize);
0924     plotinfo.legend->SetFillStyle(0);
0925     plotinfo.hstack = &hstack;
0926     plotinfo.h = plotinfo.h1 = plotinfo.h2 = 0;
0927     plotinfo.firsthisto = true;
0928 
0929     openSummaryFile();
0930     vmean.clear();
0931     vrms.clear();
0932     vdeltamean.clear();
0933     vmeanerror.clear();
0934     vPValueEqualSplitMeans.clear(), vAlignmentUncertainty.clear();
0935     vPValueMeanEqualIdeal.clear();
0936     vPValueRMSEqualIdeal.clear();
0937 
0938     std::string stringsubdet;
0939     switch (i) {
0940       case 1:
0941         stringsubdet = "BPIX";
0942         break;
0943       case 2:
0944         stringsubdet = "FPIX";
0945         break;
0946       case 3:
0947         stringsubdet = "TIB";
0948         break;
0949       case 4:
0950         stringsubdet = "TID";
0951         break;
0952       case 5:
0953         stringsubdet = "TOB";
0954         break;
0955       case 6:
0956         stringsubdet = "TEC";
0957         break;
0958     }
0959 
0960     for (std::vector<TkOfflineVariables*>::iterator it = sourceList.begin(); it != sourceList.end(); ++it) {
0961       plotinfo.vars = *it;
0962       plotinfo.h1 = plotinfo.h2 = plotinfo.h = 0;
0963 
0964       int minlayer = plotLayers ? 1 : plotLayerN;
0965       //Layer 0 is associated to the entire structure, this check ensures that even when both the plotLayers and the plotPlain options are active, also the histogram for the entire structure is made.
0966       if (plotinfo.plotPlain)
0967         minlayer = 0;
0968       int maxlayer = plotLayers ? numberOfLayers(plotinfo.vars->getPhase(), plotinfo.subDetId) : plotLayerN;
0969 
0970       for (int layer = minlayer; layer <= maxlayer; layer++) {
0971         if (plotinfo.plotPlain) {
0972           plotDMRHistogram(plotinfo, 0, layer, stringsubdet);
0973         }
0974 
0975         if (plotinfo.plotSplits) {
0976           plotDMRHistogram(plotinfo, -1, layer, stringsubdet);
0977           plotDMRHistogram(plotinfo, 1, layer, stringsubdet);
0978         }
0979 
0980         if (plotinfo.plotPlain) {
0981           if (plotinfo.h) {
0982             setDMRHistStyleAndLegend(plotinfo.h, plotinfo, 0, layer);
0983           } else {
0984             if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") &&
0985                 /*!plotinfo.plotLayers && */ layer == 0) {
0986               vmean.push_back(nan(""));
0987               vrms.push_back(nan(""));
0988               vmeanerror.push_back(nan(""));
0989               vAlignmentUncertainty.push_back(nan(""));
0990               vPValueMeanEqualIdeal.push_back(nan(""));
0991               vPValueRMSEqualIdeal.push_back(nan(""));
0992               if (plotinfo.plotSplits && plotinfo.plotPlain) {
0993                 vdeltamean.push_back(nan(""));
0994                 vPValueEqualSplitMeans.push_back(nan(""));
0995               }
0996             }
0997           }
0998         }
0999 
1000         if (plotinfo.plotSplits) {
1001           // Add delta mu to the histogram
1002           if (plotinfo.h1 != 0 && plotinfo.h2 != 0 && !plotinfo.plotPlain) {
1003             std::ostringstream legend;
1004             std::string unit = " #mum";
1005             legend.precision(3);
1006             legend << fixed;  // to always show 3 decimals
1007             float factor = 10000.0f;
1008             if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" ||
1009                 plotinfo.variable == "rmsNormX" || plotinfo.variable == "rmsNormY") {
1010               factor = 1.0f;
1011               unit = "";
1012             }
1013             float deltamu = factor * (plotinfo.h2->GetMean(1) - plotinfo.h1->GetMean(1));
1014             legend << plotinfo.vars->getName();
1015             if (layer > 0) {
1016               // TEC and TID have discs, the rest have layers
1017               if (i == 4 || i == 6)
1018                 legend << ", disc ";
1019               else
1020                 legend << ", layer ";
1021               legend << layer;
1022             }
1023             plotinfo.legend->AddEntry(static_cast<TObject*>(0), legend.str().c_str(), "");
1024             legend.str("");
1025             legend << "#Delta#mu = " << deltamu << unit;
1026             plotinfo.legend->AddEntry(static_cast<TObject*>(0), legend.str().c_str(), "");
1027 
1028             if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && !plotLayers && layer == 0) {
1029               vdeltamean.push_back(deltamu);
1030             }
1031           }
1032           if (plotinfo.h1) {
1033             setDMRHistStyleAndLegend(plotinfo.h1, plotinfo, -1, layer);
1034           }
1035           if (plotinfo.h2) {
1036             setDMRHistStyleAndLegend(plotinfo.h2, plotinfo, 1, layer);
1037           }
1038         }
1039       }
1040     }
1041 
1042     if (hstack.GetHists() != 0 && hstack.GetHists()->GetSize() != 0) {
1043       hstack.Draw("nostack");
1044       hstack.SetMaximum(plotinfo.maxY * 1.3);
1045       setTitleStyle(hstack, variable.c_str(), "#modules", plotinfo.subDetId);
1046       setHistStyle(*hstack.GetHistogram(), variable.c_str(), "#modules", 1);
1047 
1048       plotinfo.legend->Draw();
1049     } else {
1050       // Draw an empty default histogram
1051       plotinfo.h = new TH1F("defhist", "Empty default histogram", plotinfo.nbins, plotinfo.min, plotinfo.max);
1052       plotinfo.h->SetMaximum(10);
1053       if (plotinfo.variable.find("Norm") == std::string::npos)
1054         scaleXaxis(plotinfo.h, 10000);
1055       setTitleStyle(*plotinfo.h, variable.c_str(), "#modules", plotinfo.subDetId);
1056       setHistStyle(*plotinfo.h, variable.c_str(), "#modules", 1);
1057       plotinfo.h->Draw();
1058     }
1059 
1060     std::ostringstream plotName;
1061     plotName << outputDir << "/D";
1062 
1063     if (variable == "medianX")
1064       plotName << "medianR_";
1065     else if (variable == "medianY")
1066       plotName << "medianYR_";
1067     else if (variable == "meanX")
1068       plotName << "meanR_";
1069     else if (variable == "meanY")
1070       plotName << "meanYR_";
1071     else if (variable == "meanNormX")
1072       plotName << "meanNR_";
1073     else if (variable == "meanNormY")
1074       plotName << "meanNYR_";
1075     else if (variable == "rmsX")
1076       plotName << "rmsR_";
1077     else if (variable == "rmsY")
1078       plotName << "rmsYR_";
1079     else if (variable == "rmsNormX")
1080       plotName << "rmsNR_";
1081     else if (variable == "rmsNormY")
1082       plotName << "rmsNYR_";
1083 
1084     TString subdet;
1085     switch (i) {
1086       case 1:
1087         subdet = "BPIX";
1088         break;
1089       case 2:
1090         subdet = "FPIX";
1091         break;
1092       case 3:
1093         subdet = "TIB";
1094         break;
1095       case 4:
1096         subdet = "TID";
1097         break;
1098       case 5:
1099         subdet = "TOB";
1100         break;
1101       case 6:
1102         subdet = "TEC";
1103         break;
1104     }
1105 
1106     plotName << subdet;
1107 
1108     if (plotPlain && !plotSplits) {
1109       plotName << "_plain";
1110     } else if (!plotPlain && plotSplits) {
1111       plotName << "_split";
1112     }
1113     if (plotLayers) {
1114       // TEC and TID have discs, the rest have layers
1115       if (i == 4 || i == 6)
1116         plotName << "_discs";
1117       else
1118         plotName << "_layers";
1119     }
1120     if (plotLayerN > 0) {
1121       // TEC and TID have discs, the rest have layers
1122       if (i == 4 || i == 6)
1123         plotName << "_disc";
1124       else
1125         plotName << "_layer";
1126       plotName << plotLayerN;
1127     }
1128 
1129     // PNG,EPS,PDF files
1130     c.Update();
1131     c.Print((plotName.str() + ".png").c_str());
1132     c.Print((plotName.str() + ".eps").c_str());
1133     c.Print((plotName.str() + ".pdf").c_str());
1134 
1135     // ROOT file
1136     TFile f((plotName.str() + ".root").c_str(), "recreate");
1137     c.Write();
1138     f.Close();
1139 
1140     // Free allocated memory.
1141     delete plotinfo.h;
1142     delete plotinfo.h1;
1143     delete plotinfo.h2;
1144 
1145     if (vmean.size()) {
1146       summaryfile << "   mu_" << subdet;
1147       if (plotinfo.variable == "medianY")
1148         summaryfile << "_y";
1149       summaryfile << " (um)\t"
1150                   << "latexname=$\\mu_\\text{" << subdet << "}";
1151       if (plotinfo.variable == "medianY")
1152         summaryfile << "^{y}";
1153       summaryfile << "$ ($\\mu$m)\t"
1154                   << "format={:.3g}\t"
1155                   << "latexformat=${:.3g}$";
1156       for (auto mu : vmean)
1157         summaryfile << "\t" << mu;
1158       summaryfile << "\n";
1159     }
1160     if (vrms.size()) {
1161       summaryfile << "sigma_" << subdet;
1162       if (plotinfo.variable == "medianY")
1163         summaryfile << "_y";
1164       summaryfile << " (um)\t"
1165                   << "latexname=$\\sigma_\\text{" << subdet << "}";
1166       if (plotinfo.variable == "medianY")
1167         summaryfile << "^{y}";
1168       summaryfile << "$ ($\\mu$m)\t"
1169                   << "format={:.3g}\t"
1170                   << "latexformat=${:.3g}$";
1171       for (auto sigma : vrms)
1172         summaryfile << "\t" << sigma;
1173       summaryfile << "\n";
1174     }
1175     if (vdeltamean.size()) {
1176       summaryfile << "  dmu_" << subdet;
1177       if (plotinfo.variable == "medianY")
1178         summaryfile << "_y";
1179       summaryfile << " (um)\t"
1180                   << "latexname=$\\Delta\\mu_\\text{" << subdet << "}";
1181       if (plotinfo.variable == "medianY")
1182         summaryfile << "^{y}";
1183       summaryfile << "$ ($\\mu$m)\t"
1184                   << "format={:.3g}\t"
1185                   << "latexformat=${:.3g}$";
1186       for (auto dmu : vdeltamean)
1187         summaryfile << "\t" << dmu;
1188       summaryfile << "\n";
1189     }
1190     if (vmeanerror.size()) {
1191       summaryfile << "  sigma_mu_" << subdet;
1192       if (plotinfo.variable == "medianY")
1193         summaryfile << "_y";
1194       summaryfile << " (um)\t"
1195                   << "latexname=$\\sigma\\mu_\\text{" << subdet << "}";
1196       if (plotinfo.variable == "medianY")
1197         summaryfile << "^{y}";
1198       summaryfile << "$ ($\\mu$m)\t"
1199                   << "format={:.3g}\t"
1200                   << "latexformat=${:.3g}$";
1201       for (auto dmu : vmeanerror)
1202         summaryfile << "\t" << dmu;
1203       summaryfile << "\n";
1204     }
1205     if (vPValueEqualSplitMeans.size()) {
1206       summaryfile << "  p_delta_mu_equal_zero_" << subdet;
1207       if (plotinfo.variable == "medianY")
1208         summaryfile << "_y";
1209       summaryfile << "\t"
1210                   << "latexname=$P(\\delta\\mu_\\text{" << subdet << "}=0)";
1211       if (plotinfo.variable == "medianY")
1212         summaryfile << "^{y}";
1213       summaryfile << "$\t"
1214                   << "format={:.3g}\t"
1215                   << "latexformat=${:.3g}$";
1216       for (auto dmu : vPValueEqualSplitMeans)
1217         summaryfile << "\t" << dmu;
1218       summaryfile << "\n";
1219     }
1220     if (vAlignmentUncertainty.size()) {
1221       summaryfile << "  alignment_uncertainty_" << subdet;
1222       if (plotinfo.variable == "medianY")
1223         summaryfile << "_y";
1224       summaryfile << " (um)\t"
1225                   << "latexname=$\\sigma_\\text{align}_\\text{" << subdet << "}";
1226       if (plotinfo.variable == "medianY")
1227         summaryfile << "^{y}";
1228       summaryfile << "$ ($\\mu$m)\t"
1229                   << "format={:.3g}\t"
1230                   << "latexformat=${:.3g}$";
1231       for (auto dmu : vAlignmentUncertainty)
1232         summaryfile << "\t" << dmu;
1233       summaryfile << "\n";
1234     }
1235     if (vPValueMeanEqualIdeal.size()) {
1236       summaryfile << "  p_mean_equals_ideal_" << subdet;
1237       if (plotinfo.variable == "medianY")
1238         summaryfile << "_y";
1239       summaryfile << "\t"
1240                   << "latexname=$P(\\mu_\\text{" << subdet << "}=\\mu_\\text{ideal})";
1241       if (plotinfo.variable == "medianY")
1242         summaryfile << "^{y}";
1243       summaryfile << "$\t"
1244                   << "format={:.3g}\t"
1245                   << "latexformat=${:.3g}$";
1246       for (auto dmu : vPValueMeanEqualIdeal)
1247         summaryfile << "\t" << dmu;
1248       summaryfile << "\n";
1249     }
1250     if (vPValueRMSEqualIdeal.size()) {
1251       summaryfile << "  p_RMS_equals_ideal_" << subdet;
1252       if (plotinfo.variable == "medianY")
1253         summaryfile << "_y";
1254       summaryfile << "\t"
1255                   << "latexname=$P(\\sigma_\\text{" << subdet << "}=\\sigma_\\text{ideal})";
1256       if (plotinfo.variable == "medianY")
1257         summaryfile << "^{y}";
1258       summaryfile << "$\t"
1259                   << "format={:.3g}\t"
1260                   << "latexformat=${:.3g}$";
1261       for (auto dmu : vPValueRMSEqualIdeal)
1262         summaryfile << "\t" << dmu;
1263       summaryfile << "\n";
1264     }
1265   }
1266 }
1267 
1268 //------------------------------------------------------------------------------
1269 void PlotAlignmentValidation::plotChi2(const char* inputFile) {
1270   // Opens the file (it should be OfflineValidation(Parallel)_result.root)
1271   // and reads and plots the norm_chi^2 and h_chi2Prob -distributions.
1272 
1273   Bool_t errorflag = kTRUE;
1274   TFile* fi1 = TFile::Open(inputFile, "read");
1275   TDirectoryFile* mta1 = NULL;
1276   TDirectoryFile* mtb1 = NULL;
1277   TCanvas* normchi = NULL;
1278   TCanvas* chiprob = NULL;
1279   if (fi1 != NULL) {
1280     mta1 = (TDirectoryFile*)fi1->Get("TrackerOfflineValidationStandalone");
1281     if (mta1 != NULL) {
1282       mtb1 = (TDirectoryFile*)mta1->Get("GlobalTrackVariables");
1283       if (mtb1 != NULL) {
1284         normchi = dynamic_cast<TCanvas*>(mtb1->Get("h_normchi2"));
1285         chiprob = dynamic_cast<TCanvas*>(mtb1->Get("h_chi2Prob"));
1286         if (normchi != NULL && chiprob != NULL) {
1287           errorflag = kFALSE;
1288         }
1289       }
1290     }
1291   }
1292   if (errorflag) {
1293     std::cout << "PlotAlignmentValidation::plotChi2: Can't find data from given file,"
1294               << " no chi^2-plots produced" << std::endl;
1295     return;
1296   }
1297 
1298   TLegend* legend = 0;
1299   for (auto primitive : *normchi->GetListOfPrimitives()) {
1300     legend = dynamic_cast<TLegend*>(primitive);
1301     if (legend)
1302       break;
1303   }
1304   if (legend) {
1305     openSummaryFile();
1306     summaryfile << "ntracks";
1307     for (auto alignment : sourceList) {
1308       summaryfile << "\t";
1309       TString title = alignment->getName();
1310       int color = alignment->getLineColor();
1311       int style = alignment->getLineStyle();
1312       bool foundit = false;
1313       for (auto entry : *legend->GetListOfPrimitives()) {
1314         TLegendEntry* legendentry = dynamic_cast<TLegendEntry*>(entry);
1315         assert(legendentry);
1316         TH1* h = dynamic_cast<TH1*>(legendentry->GetObject());
1317         if (!h)
1318           continue;
1319         if (legendentry->GetLabel() == title && h->GetLineColor() == color && h->GetLineStyle() == style) {
1320           foundit = true;
1321           summaryfile << h->GetEntries();
1322           break;
1323         }
1324       }
1325       if (!foundit) {
1326         summaryfile << 0;
1327       }
1328     }
1329     summaryfile << "\n";
1330   }
1331 
1332   chiprob->Draw();
1333   normchi->Draw();
1334 
1335   // PNG,EPS,PDF files
1336   normchi->Print((outputDir + "/h_normchi2.png").c_str());
1337   chiprob->Print((outputDir + "/h_chi2Prob.png").c_str());
1338   normchi->Print((outputDir + "/h_normchi2.eps").c_str());
1339   chiprob->Print((outputDir + "/h_chi2Prob.eps").c_str());
1340   normchi->Print((outputDir + "/h_normchi2.pdf").c_str());
1341   chiprob->Print((outputDir + "/h_chi2Prob.pdf").c_str());
1342 
1343   // ROOT files
1344   TFile fi2((outputDir + "/h_normchi2.root").c_str(), "recreate");
1345   normchi->Write();
1346   fi2.Close();
1347 
1348   TFile fi3((outputDir + "/h_chi2Prob.root").c_str(), "recreate");
1349   chiprob->Write();
1350   fi3.Close();
1351 
1352   delete fi1;
1353 }
1354 
1355 //------------------------------------------------------------------------------
1356 THStack* PlotAlignmentValidation::addHists(
1357     const TString& selection, const TString& residType, TLegend** myLegend, bool printModuleIds, bool validforphase0) {
1358   enum ResidType {
1359     xPrimeRes,
1360     yPrimeRes,
1361     xPrimeNormRes,
1362     yPrimeNormRes,
1363     xRes,
1364     yRes,
1365     xNormRes, /*yResNorm*/
1366     ResXvsXProfile,
1367     ResXvsYProfile,
1368     ResYvsXProfile,
1369     ResYvsYProfile
1370   };
1371   ResidType rType = xPrimeRes;
1372   if (residType == "xPrime")
1373     rType = xPrimeRes;
1374   else if (residType == "yPrime")
1375     rType = yPrimeRes;
1376   else if (residType == "xPrimeNorm")
1377     rType = xPrimeNormRes;
1378   else if (residType == "yPrimeNorm")
1379     rType = yPrimeNormRes;
1380   else if (residType == "x")
1381     rType = xRes;
1382   else if (residType == "y")
1383     rType = yRes;
1384   else if (residType == "xNorm")
1385     rType = xNormRes;
1386   // else if (residType == "yNorm") rType = yResNorm;
1387   else if (residType == "ResXvsXProfile")
1388     rType = ResXvsXProfile;
1389   else if (residType == "ResYvsXProfile")
1390     rType = ResYvsXProfile;
1391   else if (residType == "ResXvsYProfile")
1392     rType = ResXvsYProfile;
1393   else if (residType == "ResYvsYProfile")
1394     rType = ResYvsYProfile;
1395   else {
1396     std::cout << "PlotAlignmentValidation::addHists: Unknown residual type " << residType << std::endl;
1397     return 0;
1398   }
1399 
1400   cout << "PlotAlignmentValidation::addHists: using selection " << selection << endl;
1401   THStack* retHistoStack = new THStack("hstack", "");
1402   if (myLegend != 0)
1403     if (*myLegend == 0) {
1404       *myLegend = new TLegend(0.17, 0.80, 0.85, 0.88);
1405     }
1406 
1407   for (std::vector<TkOfflineVariables*>::iterator itSourceFile = sourceList.begin(); itSourceFile != sourceList.end();
1408        ++itSourceFile) {
1409     std::vector<TString> histnames;
1410 
1411     TFile* f = (*itSourceFile)->getFile();
1412     TTree* tree = (*itSourceFile)->getTree();
1413     int myLineColor = (*itSourceFile)->getLineColor();
1414     int myLineStyle = (*itSourceFile)->getLineStyle();
1415     TString myLegendName = (*itSourceFile)->getName();
1416     TH1* h = 0;         // becomes result
1417     UInt_t nEmpty = 0;  // selected, but empty hists
1418     Long64_t nentries = tree->GetEntriesFast();
1419     if (!f || !tree) {
1420       std::cout << "PlotAlignmentValidation::addHists: no tree or no file" << std::endl;
1421       return 0;
1422     }
1423 
1424     bool histnamesfilled = false;
1425     int phase = (bool)(f->Get("TrackerOfflineValidationStandalone/Pixel/P1PXBBarrel_1"));
1426     if (residType.Contains("Res") && residType.Contains("Profile")) {
1427       TString basename = TString(residType)
1428                              .ReplaceAll("Res", "p_res")
1429                              .ReplaceAll("vs", "")
1430                              .ReplaceAll("Profile", "_");  //gives e.g.: p_resXX_
1431       if (selection == "subDetId==1") {
1432         if (phase == 1)
1433           histnames.push_back(TString(basename) += "P1PXBBarrel_1");
1434         else
1435           histnames.push_back(TString(basename) += "TPBBarrel_1");
1436         histnamesfilled = true;
1437       } else if (selection == "subDetId==2") {
1438         if (phase == 1) {
1439           histnames.push_back(TString(basename) += "P1PXECEndcap_2");
1440           histnames.push_back(TString(basename) += "P1PXECEndcap_3");
1441         } else {
1442           histnames.push_back(TString(basename) += "TPEEndcap_2");
1443           histnames.push_back(TString(basename) += "TPEEndcap_3");
1444         }
1445         histnamesfilled = true;
1446       } else if (selection == "subDetId==3") {
1447         histnames.push_back(TString(basename) += "TIBBarrel_1");
1448         histnamesfilled = true;
1449       } else if (selection == "subDetId==4") {
1450         histnames.push_back(TString(basename) += "TIDEndcap_2");
1451         histnames.push_back(TString(basename) += "TIDEndcap_3");
1452         histnamesfilled = true;
1453       } else if (selection == "subDetId==5") {
1454         histnames.push_back(TString(basename) += "TOBBarrel_4");
1455         histnamesfilled = true;
1456       } else if (selection == "subDetId==6") {  //whole TEC - doesn't happen by default but easy enough to account for
1457         histnames.push_back(TString(basename) += "TECEndcap_5");
1458         histnames.push_back(TString(basename) += "TECEndcap_6");
1459         histnamesfilled = true;
1460       } else if (selection == "subDetId==6 && ring <= 4") {
1461         //There are multiple with the same name and all are needed, so give the full path.  For these TFile::Get is used later instead of FindKeyAny.
1462         for (int iEndcap = 5; iEndcap <= 6; iEndcap++)
1463           for (int iDisk = 1; iDisk <= 9; iDisk++)
1464             for (int iSide = 1; iSide <= 2; iSide++)
1465               for (int iPetal = 1; iPetal <= 8; iPetal++)
1466                 for (int iRing = 1; iRing <= 4 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9); iRing++)
1467                 //in the higher disks, the inner rings go away.  But the numbering in the file structure removes the higher numbers
1468                 // so the numbers there do not correspond to the actual ring numbers
1469                 {
1470                   stringstream s;
1471                   s << "TrackerOfflineValidationStandalone/Strip/TECEndcap_" << iEndcap << "/TECDisk_" << iDisk
1472                     << "/TECSide_" << iSide << "/TECPetal_" << iPetal << "/" << basename << "TECRing_" << iRing;
1473                   histnames.push_back(TString(s.str()));
1474                 }
1475         histnamesfilled = true;
1476       } else if (selection == "subDetId==6 && ring > 4") {
1477         //There are multiple with the same name and all are needed, so give the full path.  For these TFile::Get is used later instead of FindKeyAny.
1478         for (int iEndcap = 5; iEndcap <= 6; iEndcap++)
1479           for (int iDisk = 1; iDisk <= 9; iDisk++)
1480             for (int iSide = 1; iSide <= 2; iSide++)
1481               for (int iPetal = 1; iPetal <= 8; iPetal++)
1482                 for (int iRing = 5 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9);
1483                      iRing <= 7 - (iDisk >= 4) - (iDisk >= 7) - (iDisk >= 9);
1484                      iRing++)
1485                 //in the higher disks, the inner rings go away.  But the numbering in the file structure removes the higher numbers
1486                 // so the numbers there do not correspond to the actual ring numbers
1487                 {
1488                   stringstream s;
1489                   s << "TrackerOfflineValidationStandalone/Strip/TECEndcap_" << iEndcap << "/TECDisk_" << iDisk
1490                     << "/TECSide_" << iSide << "/TECPetal_" << iPetal << "/" << basename << "TECRing_" << iRing;
1491                   histnames.push_back(TString(s.str()));
1492                 }
1493         histnamesfilled = true;
1494       }
1495     }
1496 
1497     Long64_t nSel = 0;
1498     if (histnamesfilled && histnames.size() > 0) {
1499       nSel = (Long64_t)histnames.size();
1500     }
1501     if (!histnamesfilled) {
1502       // first loop on tree to find out which entries (i.e. modules) fulfill the selection
1503       // 'Entry$' gives the entry number in the tree
1504       nSel = tree->Draw("Entry$", selection, "goff");
1505       if (nSel == -1)
1506         return 0;  // error in selection
1507       if (nSel == 0) {
1508         std::cout << "PlotAlignmentValidation::addHists: no selected module." << std::endl;
1509         return 0;
1510       }
1511       // copy entry numbers that fulfil the selection
1512       const std::vector<double> selected(tree->GetV1(), tree->GetV1() + nSel);
1513 
1514       std::vector<double>::const_iterator iterEnt = selected.begin();
1515 
1516       // second loop on tree:
1517       // for each selected entry get the hist from the file and merge
1518       TkOffTreeVariables* treeMem = 0;  // ROOT will initialise
1519       tree->SetBranchAddress("TkOffTreeVariables", &treeMem);
1520       for (Long64_t i = 0; i < nentries; i++) {
1521         if (i < *iterEnt - 0.1               // smaller index (with tolerance): skip
1522             || iterEnt == selected.end()) {  // at the end: skip
1523           continue;
1524         } else if (TMath::Abs(i - *iterEnt) < 0.11) {
1525           ++iterEnt;  // take this entry!
1526         } else
1527           std::cout << "Must not happen: " << i << " " << *iterEnt << std::endl;
1528 
1529         tree->GetEntry(i);
1530         if (printModuleIds) {
1531           std::cout << treeMem->moduleId << ": " << treeMem->entries << " entries" << std::endl;
1532         }
1533         if (treeMem->entries <= 0) {  // little speed up: skip empty hists
1534           ++nEmpty;
1535           continue;
1536         }
1537         TString hName;
1538         switch (rType) {
1539           case xPrimeRes:
1540             hName = treeMem->histNameX.c_str();
1541             break;
1542           case yPrimeRes:
1543             hName = treeMem->histNameY.c_str();
1544             break;
1545           case xPrimeNormRes:
1546             hName = treeMem->histNameNormX.c_str();
1547             break;
1548           case yPrimeNormRes:
1549             hName = treeMem->histNameNormY.c_str();
1550             break;
1551           case xRes:
1552             hName = treeMem->histNameLocalX.c_str();
1553             break;
1554           case yRes:
1555             hName = treeMem->histNameLocalY.c_str();
1556             break;
1557           case xNormRes:
1558             hName = treeMem->histNameNormLocalX.c_str();
1559             break;
1560             /*case yResNorm:      hName = treeMem->histNameNormLocalY.c_str(); break;*/
1561           case ResXvsXProfile:
1562             hName = treeMem->profileNameResXvsX.c_str();
1563             break;
1564           case ResXvsYProfile:
1565             hName = treeMem->profileNameResXvsY.c_str();
1566             break;
1567           case ResYvsXProfile:
1568             hName = treeMem->profileNameResYvsX.c_str();
1569             break;
1570           case ResYvsYProfile:
1571             hName = treeMem->profileNameResYvsY.c_str();
1572             break;
1573         }
1574         histnames.push_back(hName);
1575       }
1576     }
1577 
1578     for (std::vector<TString>::iterator ithistname = histnames.begin(); ithistname != histnames.end(); ++ithistname) {
1579       if (phase == 0 && !validforphase0)
1580         break;
1581       TH1* newHist;
1582       if (ithistname->Contains("/")) {
1583         newHist = (TH1*)f->Get(*ithistname);
1584       } else {
1585         TKey* histKey = f->FindKeyAny(*ithistname);
1586         newHist = (histKey ? static_cast<TH1*>(histKey->ReadObj()) : 0);
1587       }
1588       if (!newHist) {
1589         std::cout << "Hist " << *ithistname << " not found in file, break loop." << std::endl;
1590         break;
1591       }
1592       if (newHist->GetEntries() == 0) {
1593         nEmpty++;
1594         continue;
1595       }
1596       newHist->SetLineColor(myLineColor);
1597       newHist->SetLineStyle(myLineStyle);
1598       if (!h) {  // first hist: clone, but rename keeping only first part of name
1599         TString name(newHist->GetName());
1600         Ssiz_t pos_ = 0;
1601         for (UInt_t i2 = 0; i2 < 3; ++i2)
1602           pos_ = name.Index("_", pos_ + 1);
1603         name = name(0, pos_);  // only up to three '_'
1604         h = static_cast<TH1*>(newHist->Clone("summed_" + name));
1605         //      TString myTitle = Form("%s: %lld modules", selection, nSel);
1606         //  h->SetTitle( myTitle );
1607       } else {  // otherwise just add
1608         h->Add(newHist);
1609       }
1610       delete newHist;
1611     }
1612 
1613     std::cout << "PlotAlignmentValidation::addHists"
1614               << "Result is merged from " << nSel - nEmpty << " hists, " << nEmpty << " hists were empty." << std::endl;
1615 
1616     if (nSel - nEmpty == 0)
1617       continue;
1618 
1619     if (myLegend != 0)
1620       (*myLegend)->AddEntry(h, myLegendName, "L");
1621 
1622     retHistoStack->Add(h);
1623   }
1624 
1625   return retHistoStack;
1626 }
1627 
1628 //------------------------------------------------------------------------------
1629 /*! \fn fitGauss
1630  *  \brief Operate a Gaussian fit to the given histogram
1631  */
1632 TF1* PlotAlignmentValidation::fitGauss(TH1* hist, int color) {
1633   //1. fits a Gauss function to the inner range of abs(2 rms)
1634   //2. repeates the Gauss fit in a 2 sigma range around mean of first fit
1635   //returns mean and sigma from fit in micron
1636   if (!hist || hist->GetEntries() < 20)
1637     return 0;
1638 
1639   float mean = hist->GetMean();
1640   float sigma = hist->GetRMS();
1641   string functionname = "gaussian_";
1642   functionname += hist->GetName();
1643   TF1* func = new TF1(functionname.c_str(), "gaus", mean - 2. * sigma, mean + 2. * sigma);
1644 
1645   func->SetLineColor(color);
1646   func->SetLineStyle(2);
1647   if (0 == hist->Fit(func, "QNR")) {  // N: do not blow up file by storing fit!
1648     mean = func->GetParameter(1);
1649     sigma = func->GetParameter(2);
1650     // second fit: three sigma of first fit around mean of first fit
1651     func->SetRange(mean - 3. * sigma, mean + 3. * sigma);
1652     // I: integral gives more correct results if binning is too wide
1653     // L: Likelihood can treat empty bins correctly (if hist not weighted...)
1654     if (0 == hist->Fit(func, "Q0ILR")) {
1655       if (hist->GetFunction(func->GetName())) {  // Take care that it is later on drawn:
1656                                                  //hist->GetFunction(func->GetName())->ResetBit(TF1::kNotDraw);
1657       }
1658     }
1659   }
1660   return func;
1661 }
1662 
1663 //------------------------------------------------------------------------------
1664 /*! \fn storeHistogramInRootfile
1665  *  \brief Store the histogram and the gaussian function resulting from the fitGauss function into a root file
1666  */
1667 void PlotAlignmentValidation::storeHistogramInRootfile(TH1* hist) {
1668   //Store histogram and fit function in the root summary file
1669   rootsummaryfile->cd();
1670   hist->Write();
1671 }
1672 
1673 //------------------------------------------------------------------------------
1674 void PlotAlignmentValidation::scaleXaxis(TH1* hist, Int_t scale) {
1675   Double_t xmin = hist->GetXaxis()->GetXmin();
1676   Double_t xmax = hist->GetXaxis()->GetXmax();
1677   hist->GetXaxis()->SetLimits(xmin * scale, xmax * scale);
1678 }
1679 
1680 //------------------------------------------------------------------------------
1681 TObject* PlotAlignmentValidation::findObjectFromCanvas(TCanvas* canv, const char* className, Int_t n) {
1682   // Finds the n-th instance of the given class from the canvas
1683   TIter next(canv->GetListOfPrimitives());
1684   TObject* obj = 0;
1685   Int_t found = 0;
1686   while ((obj = next())) {
1687     if (strncmp(obj->ClassName(), className, 10) == 0) {
1688       if (++found == n)
1689         return obj;
1690     }
1691   }
1692 
1693   return 0;
1694 }
1695 
1696 //------------------------------------------------------------------------------
1697 void PlotAlignmentValidation::setTitleStyle(
1698     TNamed& hist, const char* titleX, const char* titleY, int subDetId, bool isSurfaceDeformation, TString secondline) {
1699   std::stringstream title_Xaxis;
1700   std::stringstream title_Yaxis;
1701   TString titleXAxis = titleX;
1702   TString titleYAxis = titleY;
1703   if (titleXAxis != "" && titleYAxis != "")
1704     cout << "plot " << titleXAxis << " vs " << titleYAxis << endl;
1705 
1706   hist.SetTitle("");
1707   TkAlStyle::drawStandardTitle();
1708 
1709   //Thanks Candice!
1710   TString subD;
1711   switch (subDetId) {
1712     case 1:
1713       subD = "BPIX";
1714       break;
1715     case 2:
1716       subD = "FPIX";
1717       break;
1718     case 3:
1719       subD = "TIB";
1720       break;
1721     case 4:
1722       subD = "TID";
1723       break;
1724     case 5:
1725       subD = "TOB";
1726       break;
1727     case 6:
1728       subD = "TEC";
1729       break;
1730   }
1731 
1732   TPaveText* text2;
1733   if (!isSurfaceDeformation) {
1734     text2 = new TPaveText(0.7, 0.3, 0.9, 0.6, "brNDC");
1735   } else {
1736     cout << "Surface Deformation" << endl;
1737     text2 = new TPaveText(0.8, 0.75, 0.9, 0.9, "brNDC");
1738   }
1739   text2->SetTextSize(0.06);
1740   text2->SetTextFont(42);
1741   text2->SetFillStyle(0);
1742   text2->SetBorderSize(0);
1743   text2->SetMargin(0.01);
1744   text2->SetTextAlign(12);  // align left
1745   text2->AddText(0.01, 0.75, subD);
1746   if (secondline != "") {
1747     text2->AddText(0.01, 0.25, secondline);
1748   }
1749   text2->Draw();
1750 }
1751 
1752 //------------------------------------------------------------------------------
1753 /*! \fn 
1754  *  \brief 
1755  */
1756 void PlotAlignmentValidation::setHistStyle(TH1& hist, const char* titleX, const char* titleY, int color) {
1757   std::stringstream title_Xaxis;
1758   std::stringstream title_Yaxis;
1759   TString titleXAxis = titleX;
1760   TString titleYAxis = titleY;
1761 
1762   if (titleXAxis.Contains("Phi"))
1763     title_Xaxis << titleX << "[rad]";
1764   else if (titleXAxis.Contains("meanX"))
1765     title_Xaxis << "#LTx'_{pred}-x'_{hit}#GT[#mum]";
1766   else if (titleXAxis.Contains("meanY"))
1767     title_Xaxis << "#LTy'_{pred}-y'_{hit}#GT[#mum]";
1768   else if (titleXAxis.Contains("rmsX"))
1769     title_Xaxis << "RMS(x'_{pred}-x'_{hit})[#mum]";
1770   else if (titleXAxis.Contains("rmsY"))
1771     title_Xaxis << "RMS(y'_{pred}-y'_{hit})[#mum]";
1772   else if (titleXAxis.Contains("meanNormX"))
1773     title_Xaxis << "#LTx'_{pred}-x'_{hit}/#sigma#GT";
1774   else if (titleXAxis.Contains("meanNormY"))
1775     title_Xaxis << "#LTy'_{pred}-y'_{hit}/#sigma#GT";
1776   else if (titleXAxis.Contains("rmsNormX"))
1777     title_Xaxis << "RMS(x'_{pred}-x'_{hit}/#sigma)";
1778   else if (titleXAxis.Contains("rmsNormY"))
1779     title_Xaxis << "RMS(y'_{pred}-y'_{hit}/#sigma)";
1780   else if (titleXAxis.Contains("meanLocalX"))
1781     title_Xaxis << "#LTx_{pred}-x_{hit}#GT[#mum]";
1782   else if (titleXAxis.Contains("rmsLocalX"))
1783     title_Xaxis << "RMS(x_{pred}-x_{hit})[#mum]";
1784   else if (titleXAxis.Contains("meanNormLocalX"))
1785     title_Xaxis << "#LTx_{pred}-x_{hit}/#sigma#GT[#mum]";
1786   else if (titleXAxis.Contains("rmsNormLocalX"))
1787     title_Xaxis << "RMS(x_{pred}-x_{hit}/#sigma)[#mum]";
1788   else if (titleXAxis.Contains("medianX"))
1789     title_Xaxis << "median(x'_{pred}-x'_{hit})[#mum]";
1790   else if (titleXAxis.Contains("medianY"))
1791     title_Xaxis << "median(y'_{pred}-y'_{hit})[#mum]";
1792   else
1793     title_Xaxis << titleX << "[cm]";
1794 
1795   if (hist.IsA()->InheritsFrom(TH1F::Class()))
1796     hist.SetLineColor(color);
1797   if (hist.IsA()->InheritsFrom(TProfile::Class())) {
1798     hist.SetMarkerStyle(20);
1799     hist.SetMarkerSize(0.8);
1800     hist.SetMarkerColor(color);
1801   }
1802 
1803   hist.GetXaxis()->SetTitle((title_Xaxis.str()).c_str());
1804 
1805   double binning = (hist.GetXaxis()->GetXmax() - hist.GetXaxis()->GetXmin()) / hist.GetNbinsX();
1806   title_Yaxis.precision(2);
1807 
1808   if (((titleYAxis.Contains("layer") || titleYAxis.Contains("ring")) && titleYAxis.Contains("subDetId")) ||
1809       titleYAxis.Contains("#modules")) {
1810     title_Yaxis << "number of modules";
1811     if (TString(title_Xaxis.str()).Contains("[#mum]"))
1812       title_Yaxis << " / " << binning << " #mum";
1813     else if (TString(title_Xaxis.str()).Contains("[cm]"))
1814       title_Yaxis << " / " << binning << " cm";
1815     else
1816       title_Yaxis << " / " << binning;
1817   } else
1818     title_Yaxis << titleY << "[cm]";
1819 
1820   hist.GetYaxis()->SetTitle((title_Yaxis.str()).c_str());
1821 
1822   hist.GetXaxis()->SetTitleFont(42);
1823   hist.GetYaxis()->SetTitleFont(42);
1824 }
1825 
1826 //------------------------------------------------------------------------------
1827 std::string PlotAlignmentValidation::getSelectionForDMRPlot(int minHits, int subDetId, int direction, int layer) {
1828   std::ostringstream builder;
1829   builder << "entries >= " << minHits;
1830   builder << " && subDetId == " << subDetId;
1831   if (direction != 0) {
1832     if (subDetId == 2) {  // FPIX is split by zDirection
1833       builder << " && zDirection == " << direction;
1834     } else {
1835       builder << " && rDirection == " << direction;
1836     }
1837   }
1838   if (layer > 0) {
1839     builder << " && layer == " << layer;
1840   }
1841   return builder.str();
1842 }
1843 
1844 std::string PlotAlignmentValidation::getVariableForDMRPlot(
1845     const std::string& histoname, const std::string& variable, int nbins, double min, double max) {
1846   std::ostringstream builder;
1847   builder << variable << ">>" << histoname << "(" << nbins << "," << min << "," << max << ")";
1848   return builder.str();
1849 }
1850 
1851 //------------------------------------------------------------------------------
1852 void PlotAlignmentValidation::setDMRHistStyleAndLegend(TH1F* h,
1853                                                        PlotAlignmentValidation::DMRPlotInfo& plotinfo,
1854                                                        int direction,
1855                                                        int layer) {
1856   TF1* fitResults = 0;
1857 
1858   h->SetDirectory(0);
1859 
1860   // The whole DMR plot is plotted with wider line than the split plots
1861   // If only split plots are plotted, they will be stronger too, though
1862   h->SetLineWidth((direction == 0 || (plotinfo.plotSplits && !plotinfo.plotPlain)) ? 2 : 1);
1863 
1864   // These lines determine the style of the plots according to rules:
1865   // -If the plot is for direction != 0, +1 or +2 is added to the given style for distinction
1866   // -However if only direction split plots are to be plotted, the additions should be 0 and +1 respectively
1867   // -Modulo 4 arithmetic, because the styles run from 1..4
1868   int linestyle = plotinfo.vars->getLineStyle() - 1, linestyleplus = 0;
1869   if (direction == 1) {
1870     linestyleplus = 1;
1871   }
1872   if (direction == -1) {
1873     linestyleplus = 2;
1874   }
1875   if (direction != 0 && plotinfo.plotSplits && !plotinfo.plotPlain) {
1876     linestyleplus--;
1877   }
1878   linestyle = (linestyle + linestyleplus) % 4 + 1;
1879 
1880   int linecolor = plotinfo.vars->getLineColor();
1881   if (plotinfo.plotLayers && layer > 0) {
1882     linecolor += layer - 1;
1883   }
1884 
1885   if (plotinfo.firsthisto) {
1886     setHistStyle(*h, plotinfo.variable.c_str(), "#modules", 1);  //set color later
1887     plotinfo.firsthisto = false;
1888   }
1889 
1890   h->SetLineColor(linecolor);
1891   h->SetLineStyle(linestyle);
1892 
1893   if (plotinfo.maxY < h->GetMaximum()) {
1894     plotinfo.maxY = h->GetMaximum();
1895   }
1896 
1897   //fit histogram for median and mean
1898   if (plotinfo.variable == "medianX" || plotinfo.variable == "meanX" || plotinfo.variable == "medianY" ||
1899       plotinfo.variable == "meanY") {
1900     fitResults = fitGauss(h, linecolor);
1901   }
1902 
1903   plotinfo.hstack->Add(h);
1904 
1905   std::ostringstream legend;
1906   legend.precision(3);
1907   legend << fixed;  // to always show 3 decimals
1908 
1909   // Legend: header part
1910   if (direction == -1 && plotinfo.subDetId != 2) {
1911     legend << "rDirection < 0";
1912   } else if (direction == 1 && plotinfo.subDetId != 2) {
1913     legend << "rDirection > 0";
1914   } else if (direction == -1 && plotinfo.subDetId == 2) {
1915     legend << "zDirection < 0";
1916   } else if (direction == 1 && plotinfo.subDetId == 2) {
1917     legend << "zDirection > 0";
1918   } else {
1919     legend << plotinfo.vars->getName();
1920     if (layer > 0) {
1921       // TEC and TID have discs, the rest have layers
1922       if (plotinfo.subDetId == 4 || plotinfo.subDetId == 6)
1923         legend << ", disc ";
1924       else
1925         legend << ", layer ";
1926       legend << layer << "";
1927     }
1928   }
1929 
1930   plotinfo.legend->AddEntry(h, legend.str().c_str(), "l");
1931   legend.str("");
1932 
1933   // Legend: Statistics
1934   double mean, meanerror, rms, rmserror;
1935   TString rmsname, units;
1936   bool showdeltamu =
1937       (plotinfo.h1 != 0 && plotinfo.h2 != 0 && plotinfo.plotSplits && plotinfo.plotPlain && direction == 0);
1938   if (plotinfo.variable == "medianX" || plotinfo.variable == "meanX" || plotinfo.variable == "medianY" ||
1939       plotinfo.variable == "meanY" || plotinfo.variable == "rmsX" || plotinfo.variable == "rmsY") {
1940     if (useFit_ && fitResults) {
1941       mean = fitResults->GetParameter(1) * 10000;
1942       meanerror = fitResults->GetParError(1) * 10000;
1943       rms = fitResults->GetParameter(2) * 10000;
1944       rmserror = fitResults->GetParError(2) * 10000;
1945       rmsname = "#sigma";
1946       delete fitResults;
1947     } else {
1948       mean = h->GetMean(1) * 10000;
1949       meanerror = h->GetMeanError(1) * 10000;
1950       rms = h->GetRMS(1) * 10000;
1951       rmserror = h->GetRMSError(1) * 10000;
1952       rmsname = "rms";
1953     }
1954     units = " #mum";
1955   } else if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" || plotinfo.variable == "rmsNormX" ||
1956              plotinfo.variable == "rmsNormY") {
1957     mean = h->GetMean(1);
1958     meanerror = h->GetMeanError(1);
1959     rms = h->GetRMS(1);
1960     rmserror = h->GetRMSError(1);
1961     rmsname = "rms";
1962     units = "";
1963   }
1964   if (showMean_) {
1965     legend << " #mu = " << mean;
1966     if (showMeanError_)
1967       legend << " #pm " << meanerror;
1968     legend << units;
1969     if (showRMS_ || showdeltamu || ((showModules_ || showUnderOverFlow_) && !twolines_))
1970       legend << ", ";
1971   }
1972   if (showRMS_) {
1973     legend << " " << rmsname << " = " << rms;
1974     if (showRMSError_)
1975       legend << " #pm " << rmserror;
1976     legend << units;
1977     if (showdeltamu || ((showModules_ || showUnderOverFlow_) && !twolines_))
1978       legend << ", ";
1979   }
1980 
1981   if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && /*!plotinfo.plotLayers && */ layer == 0 &&
1982       direction == 0) {
1983     vmean.push_back(mean);
1984     vrms.push_back(rms);
1985     vmeanerror.push_back(meanerror);
1986     TH1F* ideal = (TH1F*)plotinfo.hstack->GetHists()->At(0);
1987     TH1F* h = plotinfo.h;
1988     if (h->GetRMS() >= ideal->GetRMS()) {
1989       vAlignmentUncertainty.push_back(sqrt(pow(h->GetRMS(), 2) - pow(ideal->GetRMS(), 2)));
1990     } else {
1991       vAlignmentUncertainty.push_back(nan(""));
1992     }
1993     float p = (float)resampleTestOfEqualMeans(ideal, h, 10000);
1994     vPValueMeanEqualIdeal.push_back(p);
1995     p = resampleTestOfEqualRMS(ideal, h, 10000);
1996     vPValueRMSEqualIdeal.push_back(p);
1997   }
1998 
1999   // Legend: Delta mu for split plots
2000   if (showdeltamu) {
2001     float factor = 10000.0f;
2002     if (plotinfo.variable == "meanNormX" || plotinfo.variable == "meanNormY" || plotinfo.variable == "rmsNormX" ||
2003         plotinfo.variable == "rmsNormY") {
2004       factor = 1.0f;
2005     }
2006     float deltamu = factor * (plotinfo.h2->GetMean(1) - plotinfo.h1->GetMean(1));
2007     legend << "#Delta#mu = " << deltamu << units;
2008     if ((showModules_ || showUnderOverFlow_) && !twolines_)
2009       legend << ", ";
2010 
2011     if ((plotinfo.variable == "medianX" || plotinfo.variable == "medianY") && /*!plotinfo.plotLayers && */ layer == 0 &&
2012         direction == 0) {
2013       vdeltamean.push_back(deltamu);
2014       if (plotinfo.h1->GetEntries() && plotinfo.h2->GetEntries()) {
2015         float p = (float)resampleTestOfEqualMeans(plotinfo.h1, plotinfo.h2, 10000);
2016         vPValueEqualSplitMeans.push_back(p);
2017       }
2018     }
2019   }
2020 
2021   if (twolines_) {
2022     plotinfo.legend->AddEntry((TObject*)0, legend.str().c_str(), "");
2023     plotinfo.legend->AddEntry((TObject*)0, "", "");
2024     legend.str("");
2025   }
2026 
2027   if (!showUnderOverFlow_ && showModules_) {
2028     legend << (int)h->GetEntries() << " modules";
2029   }
2030   if (showUnderOverFlow_) {
2031     if (showModules_) {
2032       legend << (int)h->GetEntries() << " modules ("
2033              << (int)h->GetBinContent(0) + (int)h->GetBinContent(h->GetNbinsX() + 1) << " outside range)";
2034     } else {
2035       legend << (int)h->GetBinContent(0) + (int)h->GetBinContent(h->GetNbinsX() + 1) << " modules outside range";
2036     }
2037   }
2038   plotinfo.legend->AddEntry((TObject*)0, legend.str().c_str(), "");
2039 
2040   // Scale the x-axis (cm to um), if needed
2041   if (plotinfo.variable.find("Norm") == std::string::npos)
2042     scaleXaxis(h, 10000);
2043 }
2044 
2045 /*!
2046  * \fn plotDMRHistogram 
2047  * \brief Create the DMR histrogram using data stored in trees and store them in the plotinfo structure.
2048  */
2049 
2050 void PlotAlignmentValidation::plotDMRHistogram(PlotAlignmentValidation::DMRPlotInfo& plotinfo,
2051                                                int direction,
2052                                                int layer,
2053                                                std::string subdet) {
2054   TH1F* h = 0;
2055   //Create a name for the histogram that summarize all relevant information: name of the geometry, variable plotted, structure, layer, and whether the modules considered point inward or outward.
2056 
2057   TString histoname = "";
2058   if (plotinfo.variable == "medianX" || plotinfo.variable == "medianY")
2059     histoname = "median";
2060   else if (plotinfo.variable == "rmsNormX" || plotinfo.variable == "rmsNormY")
2061     histoname = "DrmsNR";
2062   histoname += "_";
2063   histoname += plotinfo.vars->getName();
2064   histoname.ReplaceAll(" ", "_");
2065   histoname += "_";
2066   histoname += subdet.c_str();
2067   if (plotinfo.variable == "medianY" || plotinfo.variable == "rmsNormY")
2068     histoname += "_y";
2069   if (layer != 0) {
2070     if (subdet == "TID" || subdet == "TEC")
2071       histoname += "_disc";
2072     else
2073       histoname += "_layer";
2074     histoname += to_string(layer);
2075   }
2076   if (direction == -1) {
2077     histoname += "_minus";
2078   } else if (direction == 1) {
2079     histoname += "_plus";
2080   } else {
2081     histoname += "";
2082   }
2083   std::string plotVariable =
2084       getVariableForDMRPlot(histoname.Data(), plotinfo.variable, plotinfo.nbins, plotinfo.min, plotinfo.max);
2085   std::string selection = getSelectionForDMRPlot(plotinfo.minHits, plotinfo.subDetId, direction, layer);
2086   plotinfo.vars->getTree()->Draw(plotVariable.c_str(), selection.c_str(), "goff");
2087   if (gDirectory)
2088     gDirectory->GetObject(histoname.Data(), h);
2089   if (h && h->GetEntries() > 0) {
2090     if (direction == -1) {
2091       plotinfo.h1 = h;
2092     } else if (direction == 1) {
2093       plotinfo.h2 = h;
2094     } else {
2095       plotinfo.h = h;
2096     }
2097   }
2098   if (plotinfo.variable == "medianX" || plotinfo.variable == "medianY" || plotinfo.variable == "rmsNormX" ||
2099       plotinfo.variable == "rmsNormY")
2100     storeHistogramInRootfile(h);
2101 }
2102 
2103 void PlotAlignmentValidation::modifySSHistAndLegend(THStack* hs, TLegend* legend) {
2104   // Add mean-y-values to the legend and scale the histograms.
2105 
2106   Double_t legendY = 0.80;
2107   bool hasheader = (TkAlStyle::legendheader != "");
2108   if (hasheader)
2109     legend->SetHeader(TkAlStyle::legendheader);
2110   legend->SetFillStyle(0);
2111   int legendsize = hs->GetHists()->GetSize() + hasheader;
2112 
2113   if (legendsize > 3)
2114     legendY -= 0.01 * (legendsize - 3);
2115   if (bigtext_) {
2116     legendY -= 0.05;
2117   }
2118   if (legendY < 0.6) {
2119     std::cerr << "Warning: Huge legend!" << std::endl;
2120     legendY = 0.6;
2121   }
2122   legend->SetY1(legendY);
2123   if (bigtext_)
2124     legend->SetTextSize(TkAlStyle::textSize);
2125 
2126   // Loop over all profiles
2127   TProfile* prof = 0;
2128   TIter next(hs->GetHists());
2129   Int_t index = hasheader;  //if hasheader, first entry is the header itself
2130   while ((prof = (TProfile*)next())) {
2131     //Scaling: from cm to um
2132     Double_t scale = 10000;
2133     prof->Scale(scale);
2134 
2135     Double_t stats[6] = {0};
2136     prof->GetStats(stats);
2137 
2138     std::ostringstream legendtext;
2139     legendtext.precision(3);
2140     legendtext << fixed;  // to always show 3 decimals
2141     legendtext << ": y mean = " << stats[4] / stats[0] * scale << " #mum";
2142 
2143     TLegendEntry* entry = (TLegendEntry*)legend->GetListOfPrimitives()->At(index);
2144     if (entry == 0)
2145       cout << "PlotAlignmentValidation::PlotAlignmentValidation::modifySSLegend: Bad legend!" << endl;
2146     else
2147       entry->SetLabel((entry->GetLabel() + legendtext.str()).c_str());
2148     index++;
2149   }
2150 
2151   // Make some room for the legend
2152   hs->SetMaximum(hs->GetMaximum("nostack PE") * 1.3);
2153 }
2154 
2155 //random variable: \sigma_{X_1}-\sigma_{X_2}-\delta_{RMS}
2156 //is centered approx around 0
2157 //null hypothesis: \delta_{RMS}=0
2158 //so \delta_\sigma is a realization of this random variable
2159 //how probable is it to get our value of \delta_\sigma?
2160 //->p-value
2161 double PlotAlignmentValidation::resampleTestOfEqualRMS(TH1F* h1, TH1F* h2, int numSamples) {
2162   //vector to store realizations of random variable
2163   vector<double> diff;
2164   diff.clear();
2165   //"true" (in bootstrap terms) difference of the samples' RMS
2166   double rmsdiff = abs(h1->GetRMS() - h2->GetRMS());
2167   //means of the samples to calculate RMS
2168   double m1 = h1->GetMean();
2169   double m2 = h2->GetMean();
2170   //realization of random variable
2171   double d1 = 0;
2172   double d2 = 0;
2173   //mean of random variable
2174   double test_mean = 0;
2175   for (int i = 0; i < numSamples; i++) {
2176     d1 = 0;
2177     d2 = 0;
2178     for (int i = 0; i < h1->GetEntries(); i++) {
2179       d1 += h1->GetRandom() - m1;
2180     }
2181     for (int i = 0; i < h2->GetEntries(); i++) {
2182       d2 += h2->GetRandom() + m2;
2183     }
2184     d1 /= h1->GetEntries();
2185     d2 /= h2->GetEntries();
2186     diff.push_back(abs(d1 - d2 - rmsdiff));
2187     test_mean += abs(d1 - d2 - rmsdiff);
2188   }
2189   test_mean /= numSamples;
2190   //p value
2191   double p = 0;
2192   for (double d : diff) {
2193     if (d > rmsdiff) {
2194       p += 1;
2195     }
2196   }
2197 
2198   p /= numSamples;
2199   return p;
2200 }
2201 
2202 //random variable: (\overline{X_1}-\mu_1)-(\overline{X_2}-\mu_2)
2203 //is centered approx around 0
2204 //null hypothesis: \mu_1-\mu_2=0
2205 //so \delta_\mu is a realization of this random variable
2206 //how probable is it to get our value of \delta_\mu?
2207 //->p-value
2208 double PlotAlignmentValidation::resampleTestOfEqualMeans(TH1F* h1, TH1F* h2, int numSamples) {
2209   //vector to store realization of random variable
2210   vector<double> diff;
2211   diff.clear();
2212   //"true" (in bootstrap terms) difference of the samples' means
2213   double meandiff = abs(h1->GetMean() - h2->GetMean());
2214   //realization of random variable
2215   double d1 = 0;
2216   double d2 = 0;
2217   //mean of random variable
2218   double test_mean = 0;
2219   for (int i = 0; i < numSamples; i++) {
2220     d1 = 0;
2221     d2 = 0;
2222     for (int i = 0; i < h1->GetEntries(); i++) {
2223       d1 += h1->GetRandom();
2224     }
2225     for (int i = 0; i < h2->GetEntries(); i++) {
2226       d2 += h2->GetRandom();
2227     }
2228     d1 /= h1->GetEntries();
2229     d2 /= h2->GetEntries();
2230     diff.push_back(abs(d1 - d2 - meandiff));
2231     test_mean += abs(d1 - d2 - meandiff);
2232   }
2233   test_mean /= numSamples;
2234   //p-value
2235   double p = 0;
2236   for (double d : diff) {
2237     if (d > meandiff) {
2238       p += 1;
2239     }
2240   }
2241 
2242   p /= numSamples;
2243   return p;
2244 }
2245 
2246 float PlotAlignmentValidation::twotailedStudentTTestEqualMean(float t, float v) {
2247   return 2 * (1 - ROOT::Math::tdistribution_cdf(abs(t), v));
2248 }
2249 
2250 const TString PlotAlignmentValidation::summaryfilename = "OfflineValidationSummary";
2251 
2252 vector<TH1*> PlotAlignmentValidation::findmodule(TFile* f, unsigned int moduleid) {
2253   //TFile *f = TFile::Open(filename, "READ");
2254   TString histnamex;
2255   TString histnamey;
2256   //read necessary branch/folder
2257   auto t = (TTree*)f->Get("TrackerOfflineValidationStandalone/TkOffVal");
2258 
2259   TkOffTreeVariables* variables = 0;
2260   t->SetBranchAddress("TkOffTreeVariables", &variables);
2261   unsigned int number_of_entries = t->GetEntries();
2262   for (unsigned int i = 0; i < number_of_entries; i++) {
2263     t->GetEntry(i);
2264     if (variables->moduleId == moduleid) {
2265       histnamex = variables->histNameX;
2266       histnamey = variables->histNameY;
2267       break;
2268     }
2269   }
2270 
2271   vector<TH1*> h;
2272 
2273   auto h1 = (TH1*)f->FindObjectAny(histnamex);
2274   auto h2 = (TH1*)f->FindObjectAny(histnamey);
2275 
2276   h1->SetDirectory(0);
2277   h2->SetDirectory(0);
2278 
2279   h.push_back(h1);
2280   h.push_back(h2);
2281 
2282   return h;
2283 }
2284 
2285 void PlotAlignmentValidation::residual_by_moduleID(unsigned int moduleid) {
2286   TCanvas* cx = new TCanvas("x_residual");
2287   TCanvas* cy = new TCanvas("y_residual");
2288   TLegend* legendx = new TLegend(0.55, 0.7, 1, 0.9);
2289   TLegend* legendy = new TLegend(0.55, 0.7, 1, 0.9);
2290 
2291   legendx->SetTextSize(0.016);
2292   legendx->SetTextAlign(12);
2293   legendy->SetTextSize(0.016);
2294   legendy->SetTextAlign(12);
2295 
2296   for (auto it : sourceList) {
2297     TFile* file = it->getFile();
2298     int color = it->getLineColor();
2299     int linestyle = it->getLineStyle();  //this you set by doing h->SetLineStyle(linestyle)
2300     TString legendname = it->getName();  //this goes in the legend
2301     vector<TH1*> hist = findmodule(file, moduleid);
2302 
2303     TString histnamex = legendname + " NEntries: " + TString(to_string(int(hist[0]->GetEntries())));
2304     hist[0]->SetTitle(histnamex);
2305     hist[0]->SetStats(0);
2306     hist[0]->Rebin(50);
2307     hist[0]->SetBit(TH1::kNoTitle);
2308     hist[0]->SetLineColor(color);
2309     hist[0]->SetLineStyle(linestyle);
2310     cx->cd();
2311     hist[0]->Draw("Same");
2312     legendx->AddEntry(hist[0], histnamex, "l");
2313 
2314     TString histnamey = legendname + " NEntries: " + TString(to_string(int(hist[1]->GetEntries())));
2315     hist[1]->SetTitle(histnamey);
2316     hist[1]->SetStats(0);
2317     hist[1]->Rebin(50);
2318     hist[1]->SetBit(TH1::kNoTitle);
2319     hist[1]->SetLineColor(color);
2320     hist[1]->SetLineStyle(linestyle);
2321     cy->cd();
2322     hist[1]->Draw("Same");
2323     legendy->AddEntry(hist[1], histnamey, "l");
2324   }
2325 
2326   TString filenamex = "x_residual_" + to_string(moduleid);
2327   TString filenamey = "y_residual_" + to_string(moduleid);
2328   cx->cd();
2329   legendx->Draw();
2330   cx->SaveAs(TString(outputDir + "/") + filenamex + ".root");
2331   cx->SaveAs(TString(outputDir + "/") + filenamex + ".pdf");
2332   cx->SaveAs(TString(outputDir + "/") + filenamex + ".png");
2333   cx->SaveAs(TString(outputDir + "/") + filenamex + ".eps");
2334 
2335   cy->cd();
2336   legendy->Draw();
2337   cy->SaveAs(TString(outputDir + "/") + filenamey + ".root");
2338   cy->SaveAs(TString(outputDir + "/") + filenamey + ".pdf");
2339   cy->SaveAs(TString(outputDir + "/") + filenamey + ".png");
2340   cy->SaveAs(TString(outputDir + "/") + filenamey + ".eps");
2341 }