OfflineValidation/src/PlotAlignmentValidation.cc

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