ApeEstimatorSummary.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    ApeEstimatorSummary
// Class:      ApeEstimatorSummary
//
//
// Original Author:  Johannes Hauk,6 2-039,+41227673512,
//         Created:  Mon Oct 11 13:44:03 CEST 2010
//         Modified by: Christian Schomakers (RWTH Aachen)
// $Id: ApeEstimatorSummary.cc,v 1.14 2012/01/26 00:06:33 hauk Exp $
//
//

// system include files
#include <memory>
#include <fstream>
#include <sstream>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/one/EDAnalyzer.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/EDMException.h"

#include "CLHEP/Matrix/SymMatrix.h"

//#include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
#include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorExtendedRcd.h"
#include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"

//...............
#include "Alignment/APEEstimation/interface/TrackerSectorStruct.h"

#include "TH1.h"
#include "TString.h"
#include "TFile.h"
#include "TDirectory.h"
#include "TTree.h"
#include "TF1.h"
#include "TMath.h"
//
// class declaration
//

class ApeEstimatorSummary : public edm::one::EDAnalyzer<> {
public:
  explicit ApeEstimatorSummary(const edm::ParameterSet&);
  ~ApeEstimatorSummary() override;

private:
  void beginJob() override;
  void analyze(const edm::Event&, const edm::EventSetup&) override;
  void endJob() override;

  void openInputFile();
  void getTrackerSectorStructs();
  void bookHists();
  std::vector<double> residualErrorBinning();
  void writeHists();

  void calculateApe();

  enum ApeWeight { wInvalid, wUnity, wEntries, wEntriesOverSigmaX2 };

  // ----------member data ---------------------------
  const edm::ParameterSet parameterSet_;

  bool firstEvent = true;

  TFile* inputFile_;

  std::map<unsigned int, TrackerSectorStruct> m_tkSector_;
  unsigned int noSectors_;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
ApeEstimatorSummary::ApeEstimatorSummary(const edm::ParameterSet& iConfig)
    : parameterSet_(iConfig), inputFile_(nullptr) {}

ApeEstimatorSummary::~ApeEstimatorSummary() {}

//
// member functions
//

void ApeEstimatorSummary::openInputFile() {
  const std::string inputFileName(parameterSet_.getParameter<std::string>("InputFile"));
  std::ifstream inputFileStream;
  // Check if baseline file exists
  inputFileStream.open(inputFileName.c_str());
  if (inputFileStream.is_open()) {
    inputFile_ = new TFile(inputFileName.c_str(), "READ");
  }
  if (inputFile_) {
    edm::LogInfo("CalculateAPE") << "Input file from loop over tracks and corresponding hits sucessfully opened";
  } else {
    edm::LogError("CalculateAPE") << "There is NO input file\n"
                                  << "...APE calculation stopped. Please check path of input file name in config:\n"
                                  << "\t" << inputFileName;
    throw edm::Exception(edm::errors::Configuration, "Bad input file name");
  }
}

void ApeEstimatorSummary::getTrackerSectorStructs() {
  // At present first of the plugins registered in TFileService needs to be the one containing the normalized residual histos per sector per error bin
  TString pluginName(inputFile_->GetListOfKeys()->At(0)->GetName());

  pluginName += "/";
  TDirectory *sectorDir(nullptr), *intervalDir(nullptr);
  bool sectorBool(true);
  unsigned int iSector(1);
  for (; sectorBool; ++iSector) {
    std::stringstream sectorName, fullSectorName;
    sectorName << "Sector_" << iSector << "/";
    fullSectorName << pluginName << sectorName.str();
    TString fullName(fullSectorName.str().c_str());
    inputFile_->cd();
    sectorDir = (TDirectory*)inputFile_->TDirectory::GetDirectory(fullName);
    if (sectorDir) {
      TrackerSectorStruct tkSector;
      sectorDir->GetObject("z_name;1", tkSector.Name);
      tkSector.name = tkSector.Name->GetTitle();
      bool intervalBool(true);
      for (unsigned int iInterval(1); intervalBool; ++iInterval) {
        std::stringstream intervalName, fullIntervalName;
        intervalName << "Interval_" << iInterval << "/";
        fullIntervalName << fullSectorName.str() << intervalName.str();
        fullName = fullIntervalName.str().c_str();
        intervalDir = (TDirectory*)inputFile_->TDirectory::GetDirectory(fullName);
        if (intervalDir) {
          intervalDir->GetObject("h_sigmaX;1", tkSector.m_binnedHists[iInterval]["sigmaX"]);
          intervalDir->GetObject("h_norResX;1", tkSector.m_binnedHists[iInterval]["norResX"]);
          intervalDir->GetObject("h_sigmaY;1", tkSector.m_binnedHists[iInterval]["sigmaY"]);
          intervalDir->GetObject("h_norResY;1", tkSector.m_binnedHists[iInterval]["norResY"]);
          if (tkSector.m_binnedHists[iInterval]["sigmaY"] && tkSector.m_binnedHists[iInterval]["norResY"]) {
            tkSector.isPixel = true;
          }
        } else {
          intervalBool = false;
          if (iSector == 1)
            edm::LogInfo("CalculateAPE") << "There are " << iInterval - 1
                                         << " intervals per sector defined in input file";
        }
      }
      TDirectory* resultsDir(nullptr);
      std::stringstream fullResultName;
      fullResultName << fullSectorName.str() << "Results/";
      fullName = fullResultName.str().c_str();
      resultsDir = (TDirectory*)inputFile_->TDirectory::GetDirectory(fullName);
      if (resultsDir) {
        resultsDir->GetObject("h_entriesX;1", tkSector.EntriesX);
        if (tkSector.isPixel)
          resultsDir->GetObject("h_entriesY;1", tkSector.EntriesY);
        TTree* rawIdTree(nullptr);
        resultsDir->GetObject("rawIdTree", rawIdTree);
        unsigned int rawId(0);
        rawIdTree->SetBranchAddress("RawId", &rawId);
        for (int entry = 0; entry < rawIdTree->GetEntries(); ++entry) {
          rawIdTree->GetEntry(entry);
          // command "hadd" adds entries in TTree, so rawId are existing as often as number of files are added
          bool alreadyAdded(false);
          for (auto const& i_rawId : tkSector.v_rawId) {
            if (rawId == i_rawId)
              alreadyAdded = true;
          }
          if (alreadyAdded)
            break;
          tkSector.v_rawId.push_back(rawId);
        }
      }
      m_tkSector_[iSector] = tkSector;
    } else {
      sectorBool = false;
      edm::LogInfo("CalculateAPE") << "There are " << iSector - 1 << " sectors defined in input file";
    }
  }
  noSectors_ = iSector + 1;
}

void ApeEstimatorSummary::bookHists() {
  const std::vector<double> v_binX(this->residualErrorBinning());
  for (auto& i_sector : m_tkSector_) {
    i_sector.second.WeightX = new TH1F("h_weightX",
                                       "relative weight w_{x}/w_{tot,x};#sigma_{x}  [#mum];w_{x}/w_{tot,x}",
                                       v_binX.size() - 1,
                                       &(v_binX[0]));
    i_sector.second.MeanX = new TH1F("h_meanX",
                                     "residual mean <r_{x}/#sigma_{r,x}>;#sigma_{x}  [#mum];<r_{x}/#sigma_{r,x}>",
                                     v_binX.size() - 1,
                                     &(v_binX[0]));
    i_sector.second.RmsX = new TH1F("h_rmsX",
                                    "residual rms RMS(r_{x}/#sigma_{r,x});#sigma_{x}  [#mum];RMS(r_{x}/#sigma_{r,x})",
                                    v_binX.size() - 1,
                                    &(v_binX[0]));
    i_sector.second.FitMeanX1 = new TH1F(
        "h_fitMeanX1", "fitted residual mean #mu_{x};#sigma_{x}  [#mum];#mu_{x}", v_binX.size() - 1, &(v_binX[0]));
    i_sector.second.ResidualWidthX1 = new TH1F("h_residualWidthX1",
                                               "residual width #Delta_{x};#sigma_{x}  [#mum];#Delta_{x}",
                                               v_binX.size() - 1,
                                               &(v_binX[0]));
    i_sector.second.CorrectionX1 = new TH1F("h_correctionX1",
                                            "correction to APE_{x};#sigma_{x}  [#mum];#Delta#sigma_{align,x}  [#mum]",
                                            v_binX.size() - 1,
                                            &(v_binX[0]));
    i_sector.second.FitMeanX2 = new TH1F(
        "h_fitMeanX2", "fitted residual mean #mu_{x};#sigma_{x}  [#mum];#mu_{x}", v_binX.size() - 1, &(v_binX[0]));
    i_sector.second.ResidualWidthX2 = new TH1F("h_residualWidthX2",
                                               "residual width #Delta_{x};#sigma_{x}  [#mum];#Delta_{x}",
                                               v_binX.size() - 1,
                                               &(v_binX[0]));
    i_sector.second.CorrectionX2 = new TH1F("h_correctionX2",
                                            "correction to APE_{x};#sigma_{x}  [#mum];#Delta#sigma_{align,x}  [#mum]",
                                            v_binX.size() - 1,
                                            &(v_binX[0]));

    if (i_sector.second.isPixel) {
      i_sector.second.WeightY = new TH1F("h_weightY",
                                         "relative weight w_{y}/w_{tot,y};#sigma_{y}  [#mum];w_{y}/w_{tot,y}",
                                         v_binX.size() - 1,
                                         &(v_binX[0]));
      i_sector.second.MeanY = new TH1F("h_meanY",
                                       "residual mean <r_{y}/#sigma_{r,y}>;#sigma_{y}  [#mum];<r_{y}/#sigma_{r,y}>",
                                       v_binX.size() - 1,
                                       &(v_binX[0]));
      i_sector.second.RmsY = new TH1F("h_rmsY",
                                      "residual rms RMS(r_{y}/#sigma_{r,y});#sigma_{y}  [#mum];RMS(r_{y}/#sigma_{r,y})",
                                      v_binX.size() - 1,
                                      &(v_binX[0]));
      i_sector.second.FitMeanY1 = new TH1F(
          "h_fitMeanY1", "fitted residual mean #mu_{y};#sigma_{y}  [#mum];#mu_{y}", v_binX.size() - 1, &(v_binX[0]));
      i_sector.second.ResidualWidthY1 = new TH1F("h_residualWidthY1",
                                                 "residual width #Delta_{y};#sigma_{y}  [#mum];#Delta_{y}",
                                                 v_binX.size() - 1,
                                                 &(v_binX[0]));
      i_sector.second.CorrectionY1 = new TH1F("h_correctionY1",
                                              "correction to APE_{y};#sigma_{y}  [#mum];#Delta#sigma_{align,y}  [#mum]",
                                              v_binX.size() - 1,
                                              &(v_binX[0]));
      i_sector.second.FitMeanY2 = new TH1F(
          "h_fitMeanY2", "fitted residual mean #mu_{y};#sigma_{y}  [#mum];#mu_{y}", v_binX.size() - 1, &(v_binX[0]));
      i_sector.second.ResidualWidthY2 = new TH1F("h_residualWidthY2",
                                                 "residual width #Delta_{y};#sigma_{y}  [#mum];#Delta_{y}",
                                                 v_binX.size() - 1,
                                                 &(v_binX[0]));
      i_sector.second.CorrectionY2 = new TH1F("h_correctionY2",
                                              "correction to APE_{y};#sigma_{y}  [#mum];#Delta#sigma_{align,y}  [#mum]",
                                              v_binX.size() - 1,
                                              &(v_binX[0]));
    }
  }
}

std::vector<double> ApeEstimatorSummary::residualErrorBinning() {
  std::vector<double> v_binX;
  TH1* EntriesX(m_tkSector_[1].EntriesX);
  for (int iBin = 1; iBin <= EntriesX->GetNbinsX() + 1; ++iBin) {
    v_binX.push_back(EntriesX->GetBinLowEdge(iBin));
  }
  return v_binX;
}

void ApeEstimatorSummary::writeHists() {
  TFile* resultsFile = new TFile(parameterSet_.getParameter<std::string>("ResultsFile").c_str(), "RECREATE");
  TDirectory* baseDir = resultsFile->mkdir("ApeEstimatorSummary");
  for (auto const& i_sector : m_tkSector_) {
    std::stringstream dirName;
    dirName << "Sector_" << i_sector.first;
    TDirectory* dir = baseDir->mkdir(dirName.str().c_str());
    dir->cd();

    i_sector.second.Name->Write();

    i_sector.second.WeightX->Write();
    i_sector.second.MeanX->Write();
    i_sector.second.RmsX->Write();
    i_sector.second.FitMeanX1->Write();
    i_sector.second.ResidualWidthX1->Write();
    i_sector.second.CorrectionX1->Write();
    i_sector.second.FitMeanX2->Write();
    i_sector.second.ResidualWidthX2->Write();
    i_sector.second.CorrectionX2->Write();

    if (i_sector.second.isPixel) {
      i_sector.second.WeightY->Write();
      i_sector.second.MeanY->Write();
      i_sector.second.RmsY->Write();
      i_sector.second.FitMeanY1->Write();
      i_sector.second.ResidualWidthY1->Write();
      i_sector.second.CorrectionY1->Write();
      i_sector.second.FitMeanY2->Write();
      i_sector.second.ResidualWidthY2->Write();
      i_sector.second.CorrectionY2->Write();
    }
  }
  resultsFile->Close();
}

void ApeEstimatorSummary::calculateApe() {
  // Set baseline or calculate APE value?
  const bool setBaseline(parameterSet_.getParameter<bool>("setBaseline"));

  // Read in baseline file for calculation of APE value (if not setting baseline)
  // Has same format as iterationFile
  const std::string baselineFileName(parameterSet_.getParameter<std::string>("BaselineFile"));
  TFile* baselineFile(nullptr);
  TTree* baselineTreeX(nullptr);
  TTree* baselineTreeY(nullptr);
  TTree* sectorNameBaselineTree(nullptr);
  if (!setBaseline) {
    std::ifstream baselineFileStream;
    // Check if baseline file exists
    baselineFileStream.open(baselineFileName.c_str());
    if (baselineFileStream.is_open()) {
      baselineFileStream.close();
      baselineFile = new TFile(baselineFileName.c_str(), "READ");
    }
    if (baselineFile) {
      edm::LogInfo("CalculateAPE") << "Baseline file for APE values sucessfully opened";
      baselineFile->GetObject("iterTreeX;1", baselineTreeX);
      baselineFile->GetObject("iterTreeY;1", baselineTreeY);
      baselineFile->GetObject("nameTree;1", sectorNameBaselineTree);
    } else {
      edm::LogWarning("CalculateAPE") << "There is NO baseline file for APE values, so normalized residual width =1 "
                                         "for ideal conditions is assumed";
    }
  }

  // Set up root file for iterations on APE value (or for setting baseline in setBaseline mode)
  const std::string iterationFileName(setBaseline ? baselineFileName
                                                  : parameterSet_.getParameter<std::string>("IterationFile"));
  // For iterations, updates are needed to not overwrite the iterations before
  TFile* iterationFile = new TFile(iterationFileName.c_str(), setBaseline ? "RECREATE" : "UPDATE");

  // Set up TTree for iterative APE values on first pass (first iteration) or read from file (further iterations)
  TTree* iterationTreeX(nullptr);
  TTree* iterationTreeY(nullptr);
  iterationFile->GetObject("iterTreeX;1", iterationTreeX);
  iterationFile->GetObject("iterTreeY;1", iterationTreeY);
  // The same for TTree containing the names of the sectors (no additional check, since always handled exactly as iterationTree)
  TTree* sectorNameTree(nullptr);
  iterationFile->GetObject("nameTree;1", sectorNameTree);

  bool firstIter(false);
  if (!iterationTreeX) {  // should be always true in setBaseline mode, since file is recreated
    firstIter = true;
    if (!setBaseline) {
      iterationTreeX = new TTree("iterTreeX", "Tree for APE x values of all iterations");
      iterationTreeY = new TTree("iterTreeY", "Tree for APE y values of all iterations");
      edm::LogInfo("CalculateAPE") << "First APE iteration (number 0.), create iteration file with TTree";
      sectorNameTree = new TTree("nameTree", "Tree with names of sectors");
    } else {
      iterationTreeX = new TTree("iterTreeX", "Tree for baseline x values of normalized residual width");
      iterationTreeY = new TTree("iterTreeY", "Tree for baseline y values of normalized residual width");
      edm::LogInfo("CalculateAPE") << "Set baseline, create baseline file with TTree";
      sectorNameTree = new TTree("nameTree", "Tree with names of sectors");
    }
  } else {
    const unsigned int iteration(iterationTreeX->GetEntries());
    edm::LogWarning("CalculateAPE") << "NOT the first APE iteration (number 0.) but the " << iteration
                                    << ". one, is this wanted or forgot to delete old iteration file with TTree?";
  }

  // Assign the information stored in the trees to arrays
  double a_apeSectorX[noSectors_];
  double a_apeSectorY[noSectors_];
  double a_baselineSectorX[noSectors_];
  double a_baselineSectorY[noSectors_];

  std::string* a_sectorName[noSectors_];
  std::string* a_sectorBaselineName[noSectors_];
  for (auto const& i_sector : m_tkSector_) {
    const unsigned int iSector(i_sector.first);
    const bool pixelSector(i_sector.second.isPixel);
    a_apeSectorX[iSector] = 99.;
    a_apeSectorY[iSector] = 99.;
    a_baselineSectorX[iSector] = -99.;
    a_baselineSectorY[iSector] = -99.;

    a_sectorName[iSector] = nullptr;
    a_sectorBaselineName[iSector] = nullptr;
    std::stringstream ss_sector, ss_sectorSuffixed;
    ss_sector << "Ape_Sector_" << iSector;
    if (!setBaseline && baselineTreeX) {
      baselineTreeX->SetBranchAddress(ss_sector.str().c_str(), &a_baselineSectorX[iSector]);
      baselineTreeX->GetEntry(0);
      if (pixelSector) {
        baselineTreeY->SetBranchAddress(ss_sector.str().c_str(), &a_baselineSectorY[iSector]);
        baselineTreeY->GetEntry(0);
      }
      sectorNameBaselineTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorBaselineName[iSector]);
      sectorNameBaselineTree->GetEntry(0);
    } else {
      // Set default ideal normalized residual width to 1
      a_baselineSectorX[iSector] = 1.;
      a_baselineSectorY[iSector] = 1.;
    }
    if (firstIter) {  // should be always true in setBaseline mode, since file is recreated
      ss_sectorSuffixed << ss_sector.str() << "/D";
      iterationTreeX->Branch(ss_sector.str().c_str(), &a_apeSectorX[iSector], ss_sectorSuffixed.str().c_str());
      if (pixelSector) {
        iterationTreeY->Branch(ss_sector.str().c_str(), &a_apeSectorY[iSector], ss_sectorSuffixed.str().c_str());
      }
      sectorNameTree->Branch(ss_sector.str().c_str(), &a_sectorName[iSector], 32000, 00);
    } else {
      iterationTreeX->SetBranchAddress(ss_sector.str().c_str(), &a_apeSectorX[iSector]);
      iterationTreeX->GetEntry(iterationTreeX->GetEntries() - 1);
      if (pixelSector) {
        iterationTreeY->SetBranchAddress(ss_sector.str().c_str(), &a_apeSectorY[iSector]);
        iterationTreeY->GetEntry(iterationTreeY->GetEntries() - 1);
      }
      sectorNameTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorName[iSector]);
      sectorNameTree->GetEntry(0);
    }
  }

  // Check whether sector definitions are identical with the ones of previous iterations and with the ones in baseline file
  for (auto& i_sector : m_tkSector_) {
    const std::string& name(i_sector.second.name);
    if (!firstIter) {
      const std::string& nameLastIter(*a_sectorName[i_sector.first]);
      if (name != nameLastIter) {
        edm::LogError("CalculateAPE") << "Inconsistent sector definition in iterationFile for sector " << i_sector.first
                                      << ",\n"
                                      << "Recent iteration has name \"" << name << "\", while previous had \""
                                      << nameLastIter << "\"\n"
                                      << "...APE calculation stopped. Please check sector definitions in config!\n";
        return;
      }
    } else {
      a_sectorName[i_sector.first] = new std::string(name);
    }
    if (!setBaseline && baselineFile) {
      const std::string& nameBaseline(*a_sectorBaselineName[i_sector.first]);
      if (name != nameBaseline) {
        edm::LogError("CalculateAPE") << "Inconsistent sector definition in baselineFile for sector " << i_sector.first
                                      << ",\n"
                                      << "Recent iteration has name \"" << name << "\", while baseline had \""
                                      << nameBaseline << "\"\n"
                                      << "...APE calculation stopped. Please check sector definitions in config!\n";
        return;
      }
    }
  }

  // Set up text file for writing out APE values per module
  std::ofstream apeOutputFile;
  if (!setBaseline) {
    const std::string apeOutputFileName(parameterSet_.getParameter<std::string>("ApeOutputFile"));
    apeOutputFile.open(apeOutputFileName.c_str());
    if (apeOutputFile.is_open()) {
      edm::LogInfo("CalculateAPE") << "Text file for writing APE values successfully opened";
    } else {
      edm::LogError("CalculateAPE") << "Text file for writing APE values NOT opened,\n"
                                    << "...APE calculation stopped. Please check path of text file name in config:\n"
                                    << "\t" << apeOutputFileName;
      return;
    }
  }

  // Loop over sectors for calculating APE
  const std::string apeWeightName(parameterSet_.getParameter<std::string>("apeWeight"));
  ApeWeight apeWeight(wInvalid);
  if (apeWeightName == "unity")
    apeWeight = wUnity;
  else if (apeWeightName == "entries")
    apeWeight = wEntries;
  else if (apeWeightName == "entriesOverSigmaX2")
    apeWeight = wEntriesOverSigmaX2;
  if (apeWeight == wInvalid) {
    edm::LogError("CalculateAPE") << "Invalid parameter 'apeWeight' in cfg file: \"" << apeWeightName
                                  << "\"\nimplemented apeWeights are \"unity\", \"entries\", \"entriesOverSigmaX2\""
                                  << "\n...APE calculation stopped.";
    return;
  }
  const double minHitsPerInterval(parameterSet_.getParameter<double>("minHitsPerInterval"));
  const double sigmaFactorFit(parameterSet_.getParameter<double>("sigmaFactorFit"));
  const double correctionScaling(parameterSet_.getParameter<double>("correctionScaling"));
  const bool smoothIteration(parameterSet_.getParameter<bool>("smoothIteration"));
  const double smoothFraction(parameterSet_.getParameter<double>("smoothFraction"));
  if (smoothFraction <= 0. || smoothFraction > 1.) {
    edm::LogError("CalculateAPE") << "Incorrect parameter in cfg file,"
                                  << "\nsmoothFraction has to be in [0,1], but is " << smoothFraction
                                  << "\n...APE calculation stopped.";
    return;
  }
  for (auto& i_sector : m_tkSector_) {
    typedef std::pair<double, double> Error2AndResidualWidth2PerBin;
    typedef std::pair<double, Error2AndResidualWidth2PerBin> WeightAndResultsPerBin;
    std::vector<WeightAndResultsPerBin> v_weightAndResultsPerBinX;
    std::vector<WeightAndResultsPerBin> v_weightAndResultsPerBinY;

    double baselineWidthX2(a_baselineSectorX[i_sector.first]);
    double baselineWidthY2(a_baselineSectorY[i_sector.first]);

    // Loop over residual error bins to calculate APE for every bin
    for (auto const& i_errBins : i_sector.second.m_binnedHists) {
      std::map<std::string, TH1*> mHists = i_errBins.second;

      double entriesX = mHists["sigmaX"]->GetEntries();
      double meanSigmaX = mHists["sigmaX"]->GetMean();

      // Fitting Parameters
      double xMin = mHists["norResX"]->GetXaxis()->GetXmin();
      double xMax = mHists["norResX"]->GetXaxis()->GetXmax();
      double integralX = mHists["norResX"]->Integral();
      double meanX = mHists["norResX"]->GetMean();
      double rmsX = mHists["norResX"]->GetRMS();
      double maximumX = mHists["norResX"]->GetBinContent(mHists["norResX"]->GetMaximumBin());

      // First Gaus Fit
      TF1 funcX_1("mygausX", "gaus", xMin, xMax);
      funcX_1.SetParameters(maximumX, meanX, rmsX);
      TString fitOpt("ILERQ");  //TString fitOpt("IMR"); ("IRLL"); ("IRQ");
      if (integralX > minHitsPerInterval) {
        if (mHists["norResX"]->Fit(&funcX_1, fitOpt)) {
          edm::LogWarning("CalculateAPE") << "Fit1 did not work: " << mHists["norResX"]->Fit(&funcX_1, fitOpt);
          continue;
        }
        LogDebug("CalculateAPE") << "FitResultX1\t" << mHists["norResX"]->Fit(&funcX_1, fitOpt) << "\n";
      }
      double meanX_1 = funcX_1.GetParameter(1);
      double sigmaX_1 = funcX_1.GetParameter(2);

      // Second gaus fit
      TF1 funcX_2("mygausX2",
                  "gaus",
                  meanX_1 - sigmaFactorFit * TMath::Abs(sigmaX_1),
                  meanX_1 + sigmaFactorFit * TMath::Abs(sigmaX_1));
      funcX_2.SetParameters(funcX_1.GetParameter(0), meanX_1, sigmaX_1);
      if (integralX > minHitsPerInterval) {
        if (mHists["norResX"]->Fit(&funcX_2, fitOpt)) {
          edm::LogWarning("CalculateAPE") << "Fit2 did not work for x : " << mHists["norResX"]->Fit(&funcX_2, fitOpt);
          continue;
        }
        LogDebug("CalculateAPE") << "FitResultX2\t" << mHists["norResX"]->Fit(&funcX_2, fitOpt) << "\n";
      }
      double meanX_2 = funcX_2.GetParameter(1);
      double sigmaX_2 = funcX_2.GetParameter(2);

      // Now the same for y coordinate
      double entriesY(0.);
      double meanSigmaY(0.);
      if (i_sector.second.isPixel) {
        entriesY = mHists["sigmaY"]->GetEntries();
        meanSigmaY = mHists["sigmaY"]->GetMean();
      }

      double meanY_1(0.);
      double sigmaY_1(0.);
      double meanY_2(0.);
      double sigmaY_2(0.);
      double meanY(0.);
      double rmsY(0.);
      if (i_sector.second.isPixel) {
        // Fitting Parameters
        double yMin = mHists["norResY"]->GetXaxis()->GetXmin();
        double yMax = mHists["norResY"]->GetXaxis()->GetXmax();
        double integralY = mHists["norResY"]->Integral();
        meanY = mHists["norResY"]->GetMean();
        rmsY = mHists["norResY"]->GetRMS();
        double maximumY = mHists["norResY"]->GetBinContent(mHists["norResY"]->GetMaximumBin());

        // First Gaus Fit
        TF1 funcY_1("mygausY", "gaus", yMin, yMax);
        funcY_1.SetParameters(maximumY, meanY, rmsY);
        if (integralY > minHitsPerInterval) {
          if (mHists["norResY"]->Fit(&funcY_1, fitOpt)) {
            edm::LogWarning("CalculateAPE") << "Fit1 did not work: " << mHists["norResY"]->Fit(&funcY_1, fitOpt);
            continue;
          }
          LogDebug("CalculateAPE") << "FitResultY1\t" << mHists["norResY"]->Fit(&funcY_1, fitOpt) << "\n";
        }
        meanY_1 = funcY_1.GetParameter(1);
        sigmaY_1 = funcY_1.GetParameter(2);

        // Second gaus fit
        TF1 funcY_2("mygausY2",
                    "gaus",
                    meanY_1 - sigmaFactorFit * TMath::Abs(sigmaY_1),
                    meanY_1 + sigmaFactorFit * TMath::Abs(sigmaY_1));
        funcY_2.SetParameters(funcY_1.GetParameter(0), meanY_1, sigmaY_1);
        if (integralY > minHitsPerInterval) {
          if (mHists["norResY"]->Fit(&funcY_2, fitOpt)) {
            edm::LogWarning("CalculateAPE") << "Fit2 did not work for y : " << mHists["norResY"]->Fit(&funcY_2, fitOpt);
            continue;
          }
          LogDebug("CalculateAPE") << "FitResultY2\t" << mHists["norResY"]->Fit(&funcY_2, fitOpt) << "\n";
        }
        meanY_2 = funcY_2.GetParameter(1);
        sigmaY_2 = funcY_2.GetParameter(2);
      }

      // Fill histograms
      double fitMeanX_1(meanX_1), fitMeanX_2(meanX_2);
      double residualWidthX_1(sigmaX_1), residualWidthX_2(sigmaX_2);
      double fitMeanY_1(meanY_1), fitMeanY_2(meanY_2);
      double residualWidthY_1(sigmaY_1), residualWidthY_2(sigmaY_2);

      double correctionX2_1(-0.0010), correctionX2_2(-0.0010);
      double correctionY2_1(-0.0010), correctionY2_2(-0.0010);
      correctionX2_1 = meanSigmaX * meanSigmaX * (residualWidthX_1 * residualWidthX_1 - baselineWidthX2);
      correctionX2_2 = meanSigmaX * meanSigmaX * (residualWidthX_2 * residualWidthX_2 - baselineWidthX2);
      if (i_sector.second.isPixel) {
        correctionY2_1 = meanSigmaY * meanSigmaY * (residualWidthY_1 * residualWidthY_1 - baselineWidthY2);
        correctionY2_2 = meanSigmaY * meanSigmaY * (residualWidthY_2 * residualWidthY_2 - baselineWidthY2);
      }

      double correctionX_1 = correctionX2_1 >= 0. ? std::sqrt(correctionX2_1) : -std::sqrt(-correctionX2_1);
      double correctionX_2 = correctionX2_2 >= 0. ? std::sqrt(correctionX2_2) : -std::sqrt(-correctionX2_2);
      double correctionY_1 = correctionY2_1 >= 0. ? std::sqrt(correctionY2_1) : -std::sqrt(-correctionY2_1);
      double correctionY_2 = correctionY2_2 >= 0. ? std::sqrt(correctionY2_2) : -std::sqrt(-correctionY2_2);
      // Meanwhile, this got very bad default values, or? (negative corrections allowed)
      if (isnan(correctionX_1))
        correctionX_1 = -0.0010;
      if (isnan(correctionX_2))
        correctionX_2 = -0.0010;
      if (isnan(correctionY_1))
        correctionY_1 = -0.0010;
      if (isnan(correctionY_2))
        correctionY_2 = -0.0010;

      if (entriesX < minHitsPerInterval) {
        meanX = 0.;
        rmsX = -0.0010;
        fitMeanX_1 = 0.;
        correctionX_1 = residualWidthX_1 = -0.0010;
        fitMeanX_2 = 0.;
        correctionX_2 = residualWidthX_2 = -0.0010;
      }

      if (entriesY < minHitsPerInterval) {
        meanY = 0.;
        rmsY = -0.0010;
        fitMeanY_1 = 0.;
        correctionY_1 = residualWidthY_1 = -0.0010;
        fitMeanY_2 = 0.;
        correctionY_2 = residualWidthY_2 = -0.0010;
      }

      i_sector.second.MeanX->SetBinContent(i_errBins.first, meanX);
      i_sector.second.RmsX->SetBinContent(i_errBins.first, rmsX);

      i_sector.second.FitMeanX1->SetBinContent(i_errBins.first, fitMeanX_1);
      i_sector.second.ResidualWidthX1->SetBinContent(i_errBins.first, residualWidthX_1);
      i_sector.second.CorrectionX1->SetBinContent(i_errBins.first, correctionX_1 * 10000.);

      i_sector.second.FitMeanX2->SetBinContent(i_errBins.first, fitMeanX_2);
      i_sector.second.ResidualWidthX2->SetBinContent(i_errBins.first, residualWidthX_2);
      i_sector.second.CorrectionX2->SetBinContent(i_errBins.first, correctionX_2 * 10000.);

      if (i_sector.second.isPixel) {
        i_sector.second.MeanY->SetBinContent(i_errBins.first, meanY);
        i_sector.second.RmsY->SetBinContent(i_errBins.first, rmsY);

        i_sector.second.FitMeanY1->SetBinContent(i_errBins.first, fitMeanY_1);
        i_sector.second.ResidualWidthY1->SetBinContent(i_errBins.first, residualWidthY_1);
        i_sector.second.CorrectionY1->SetBinContent(i_errBins.first, correctionY_1 * 10000.);

        i_sector.second.FitMeanY2->SetBinContent(i_errBins.first, fitMeanY_2);
        i_sector.second.ResidualWidthY2->SetBinContent(i_errBins.first, residualWidthY_2);
        i_sector.second.CorrectionY2->SetBinContent(i_errBins.first, correctionY_2 * 10000.);
      }

      // Use result for bin only when entries>=minHitsPerInterval
      if (entriesX < minHitsPerInterval && entriesY < minHitsPerInterval)
        continue;

      double weightX(0.);
      double weightY(0.);
      if (apeWeight == wUnity) {
        weightX = 1.;
        weightY = 1.;
      } else if (apeWeight == wEntries) {
        weightX = entriesX;
        weightY = entriesY;
      } else if (apeWeight == wEntriesOverSigmaX2) {
        weightX = entriesX / (meanSigmaX * meanSigmaX);
        weightY = entriesY / (meanSigmaY * meanSigmaY);
      }

      const Error2AndResidualWidth2PerBin error2AndResidualWidth2PerBinX(meanSigmaX * meanSigmaX,
                                                                         residualWidthX_1 * residualWidthX_1);
      const WeightAndResultsPerBin weightAndResultsPerBinX(weightX, error2AndResidualWidth2PerBinX);
      if (!(entriesX < minHitsPerInterval)) {
        //Fill absolute weights
        i_sector.second.WeightX->SetBinContent(i_errBins.first, weightX);
        v_weightAndResultsPerBinX.push_back(weightAndResultsPerBinX);
      }

      const Error2AndResidualWidth2PerBin error2AndResidualWidth2PerBinY(meanSigmaY * meanSigmaY,
                                                                         residualWidthY_1 * residualWidthY_1);
      const WeightAndResultsPerBin weightAndResultsPerBinY(weightY, error2AndResidualWidth2PerBinY);
      if (!(entriesY < minHitsPerInterval)) {
        //Fill absolute weights
        i_sector.second.WeightY->SetBinContent(i_errBins.first, weightY);
        v_weightAndResultsPerBinY.push_back(weightAndResultsPerBinY);
      }
    }

    // Do the final calculations

    if (v_weightAndResultsPerBinX.empty()) {
      edm::LogError("CalculateAPE") << "NO error interval of sector " << i_sector.first
                                    << " has a valid x APE calculated,\n...so cannot set APE";
      continue;
    }

    if (i_sector.second.isPixel && v_weightAndResultsPerBinY.empty()) {
      edm::LogError("CalculateAPE") << "NO error interval of sector " << i_sector.first
                                    << " has a valid y APE calculated,\n...so cannot set APE";
      continue;
    }

    double correctionX2(999.);
    double correctionY2(999.);

    // Get sum of all weights
    double weightSumX(0.);
    for (auto const& i_apeBin : v_weightAndResultsPerBinX) {
      weightSumX += i_apeBin.first;
    }
    i_sector.second.WeightX->Scale(1 / weightSumX);
    double weightSumY(0.);
    if (i_sector.second.isPixel) {
      for (auto const& i_apeBin : v_weightAndResultsPerBinY) {
        weightSumY += i_apeBin.first;
      }
      i_sector.second.WeightY->Scale(1 / weightSumY);
    }

    if (!setBaseline) {
      // Calculate weighted mean
      bool firstIntervalX(true);
      for (auto const& i_apeBin : v_weightAndResultsPerBinX) {
        if (firstIntervalX) {
          correctionX2 = i_apeBin.first * i_apeBin.second.first * (i_apeBin.second.second - baselineWidthX2);
          firstIntervalX = false;
        } else {
          correctionX2 += i_apeBin.first * i_apeBin.second.first * (i_apeBin.second.second - baselineWidthX2);
        }
      }
      correctionX2 = correctionX2 / weightSumX;
    } else {
      double numeratorX(0.), denominatorX(0.);
      for (auto const& i_apeBin : v_weightAndResultsPerBinX) {
        numeratorX += i_apeBin.first * i_apeBin.second.first * i_apeBin.second.second;
        denominatorX += i_apeBin.first * i_apeBin.second.first;
      }
      correctionX2 = numeratorX / denominatorX;
    }

    if (i_sector.second.isPixel) {
      if (!setBaseline) {
        // Calculate weighted mean
        bool firstIntervalY(true);
        for (auto const& i_apeBin : v_weightAndResultsPerBinY) {
          if (firstIntervalY) {
            correctionY2 = i_apeBin.first * i_apeBin.second.first * (i_apeBin.second.second - baselineWidthY2);
            firstIntervalY = false;
          } else {
            correctionY2 += i_apeBin.first * i_apeBin.second.first * (i_apeBin.second.second - baselineWidthY2);
          }
        }
        correctionY2 = correctionY2 / weightSumY;
      } else {
        double numeratorY(0.), denominatorY(0.);
        for (auto const& i_apeBin : v_weightAndResultsPerBinY) {
          numeratorY += i_apeBin.first * i_apeBin.second.first * i_apeBin.second.second;
          denominatorY += i_apeBin.first * i_apeBin.second.first;
        }
        correctionY2 = numeratorY / denominatorY;
      }
    }

    if (!setBaseline) {
      // Calculate updated squared APE of current iteration
      double apeX2(999.);
      double apeY2(999.);

      // old APE value from last iteration
      if (firstIter) {
        apeX2 = 0.;
        apeY2 = 0.;
      } else {
        apeX2 = a_apeSectorX[i_sector.first];
        apeY2 = a_apeSectorY[i_sector.first];
      }
      const double apeX2old = apeX2;
      const double apeY2old = apeY2;

      // scale APE Correction with value given in cfg (not if smoothed iteration is used)
      edm::LogInfo("CalculateAPE") << "Unscaled correction x for sector " << i_sector.first << " is "
                                   << (correctionX2 > 0. ? +1. : -1.) * std::sqrt(std::fabs(correctionX2));
      if (!smoothIteration || firstIter)
        correctionX2 = correctionX2 * correctionScaling * correctionScaling;
      if (i_sector.second.isPixel) {
        edm::LogInfo("CalculateAPE") << "Unscaled correction y for sector " << i_sector.first << " is "
                                     << (correctionY2 > 0. ? +1. : -1.) * std::sqrt(std::fabs(correctionY2));
        if (!smoothIteration || firstIter)
          correctionY2 = correctionY2 * correctionScaling * correctionScaling;
      }

      // new APE value
      // smooth iteration or not?
      if (apeX2 + correctionX2 < 0.)
        correctionX2 = -apeX2;
      if (apeY2 + correctionY2 < 0.)
        correctionY2 = -apeY2;
      const double apeX2new(apeX2old + correctionX2);
      const double apeY2new(apeY2old + correctionY2);
      if (!smoothIteration || firstIter) {
        apeX2 = apeX2new;
        apeY2 = apeY2new;
      } else {
        const double apeXtmp = smoothFraction * std::sqrt(apeX2new) + (1 - smoothFraction) * std::sqrt(apeX2old);
        const double apeYtmp = smoothFraction * std::sqrt(apeY2new) + (1 - smoothFraction) * std::sqrt(apeY2old);
        apeX2 = apeXtmp * apeXtmp;
        apeY2 = apeYtmp * apeYtmp;
      }
      if (apeX2 < 0. || apeY2 < 0.)
        edm::LogError("CalculateAPE") << "\n\n\tBad APE, but why???\n\n\n";
      a_apeSectorX[i_sector.first] = apeX2;
      a_apeSectorY[i_sector.first] = apeY2;

      // Set the calculated APE spherical for all modules of strip sectors
      const double apeX(std::sqrt(apeX2));
      const double apeY(std::sqrt(apeY2));
      const double apeZ(std::sqrt(0.5 * (apeX2 + apeY2)));
      for (auto const& i_rawId : i_sector.second.v_rawId) {
        if (i_sector.second.isPixel) {
          apeOutputFile << i_rawId << " " << std::fixed << std::setprecision(5) << apeX << " " << apeY << " " << apeZ
                        << "\n";
        } else {
          apeOutputFile << i_rawId << " " << std::fixed << std::setprecision(5) << apeX << " " << apeX << " " << apeX
                        << "\n";
        }
      }
    } else {  // In setBaseline mode, just fill estimated mean value of residual width
      a_apeSectorX[i_sector.first] = correctionX2;
      a_apeSectorY[i_sector.first] = correctionY2;
    }
  }
  if (!setBaseline)
    apeOutputFile.close();

  iterationTreeX->Fill();
  iterationTreeX->Write("iterTreeX",
                        TObject::kOverwrite);  // TObject::kOverwrite needed to not produce another iterTreeX;2
  iterationTreeY->Fill();
  iterationTreeY->Write("iterTreeY",
                        TObject::kOverwrite);  // TObject::kOverwrite needed to not produce another iterTreeY;2
  if (firstIter) {
    sectorNameTree->Fill();
    sectorNameTree->Write("nameTree");
  }
  iterationFile->Close();
  delete iterationFile;

  if (baselineFile) {
    baselineFile->Close();
    delete baselineFile;
  }
}

// ------------ method called to for each event  ------------
void ApeEstimatorSummary::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  // Load APEs from the GT and write them to root files similar to the ones from calculateAPE()

  if (firstEvent) {
    // Set baseline or calculate APE value?
    const bool setBaseline(parameterSet_.getParameter<bool>("setBaseline"));

    edm::ESHandle<AlignmentErrorsExtended> alignmentErrors;
    iSetup.get<TrackerAlignmentErrorExtendedRcd>().get(alignmentErrors);

    // Read in baseline file for calculation of APE value (if not setting baseline)
    // Has same format as iterationFile
    const std::string baselineFileName(parameterSet_.getParameter<std::string>("BaselineFile"));
    TFile* baselineFile(nullptr);
    TTree* sectorNameBaselineTree(nullptr);
    if (!setBaseline) {
      std::ifstream baselineFileStream;
      // Check if baseline file exists
      baselineFileStream.open(baselineFileName.c_str());
      if (baselineFileStream.is_open()) {
        baselineFileStream.close();
        baselineFile = new TFile(baselineFileName.c_str(), "READ");
      }
      if (baselineFile) {
        edm::LogInfo("CalculateAPE") << "Baseline file for APE values sucessfully opened";
        baselineFile->GetObject("nameTree;1", sectorNameBaselineTree);
      } else {
        edm::LogWarning("CalculateAPE") << "There is NO baseline file for APE values, so normalized residual width =1 "
                                           "for ideal conditions is assumed";
      }
    }

    // Set up root file for default APE values
    const std::string defaultFileName(parameterSet_.getParameter<std::string>("DefaultFile"));
    TFile* defaultFile = new TFile(defaultFileName.c_str(), "RECREATE");

    // Naming in the root files has to be iterTreeX to be consistent for the plotting tool
    TTree* defaultTreeX(nullptr);
    TTree* defaultTreeY(nullptr);
    defaultFile->GetObject("iterTreeX;1", defaultTreeX);
    defaultFile->GetObject("iterTreeY;1", defaultTreeY);
    // The same for TTree containing the names of the sectors (no additional check, since always handled exactly as defaultTree)
    TTree* sectorNameTree(nullptr);
    defaultFile->GetObject("nameTree;1", sectorNameTree);

    bool firstIter(false);
    if (!defaultTreeX) {  // should be always true in setBaseline mode, since file is recreated
      firstIter = true;
      defaultTreeX = new TTree("iterTreeX", "Tree for default APE x values from GT");
      defaultTreeY = new TTree("iterTreeY", "Tree for default APE y values from GT");
      edm::LogInfo("CalculateAPE") << "First APE iteration (number 0.), create default file with TTree";
      sectorNameTree = new TTree("nameTree", "Tree with names of sectors");
    } else {
      edm::LogWarning("CalculateAPE") << "NOT the first APE iteration (number 0.), is this wanted or forgot to delete "
                                         "old iteration file with TTree?";
    }

    // Assign the information stored in the trees to arrays
    double a_defaultSectorX[noSectors_];
    double a_defaultSectorY[noSectors_];

    std::string* a_sectorName[noSectors_];
    std::string* a_sectorBaselineName[noSectors_];
    for (auto const& i_sector : m_tkSector_) {
      const unsigned int iSector(i_sector.first);
      const bool pixelSector(i_sector.second.isPixel);

      a_defaultSectorX[iSector] = -99.;
      a_defaultSectorY[iSector] = -99.;

      a_sectorName[iSector] = nullptr;
      a_sectorBaselineName[iSector] = nullptr;
      std::stringstream ss_sector, ss_sectorSuffixed;
      ss_sector << "Ape_Sector_" << iSector;
      if (!setBaseline && sectorNameBaselineTree) {
        sectorNameBaselineTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorBaselineName[iSector]);
        sectorNameBaselineTree->GetEntry(0);
      }

      if (firstIter) {  // should be always true in setBaseline mode, since file is recreated
        ss_sectorSuffixed << ss_sector.str() << "/D";
        defaultTreeX->Branch(ss_sector.str().c_str(), &a_defaultSectorX[iSector], ss_sectorSuffixed.str().c_str());
        if (pixelSector) {
          defaultTreeY->Branch(ss_sector.str().c_str(), &a_defaultSectorY[iSector], ss_sectorSuffixed.str().c_str());
        }
        sectorNameTree->Branch(ss_sector.str().c_str(), &a_sectorName[iSector], 32000, 00);
      } else {
        defaultTreeX->SetBranchAddress(ss_sector.str().c_str(), &a_defaultSectorX[iSector]);
        defaultTreeX->GetEntry(0);
        if (pixelSector) {
          defaultTreeY->SetBranchAddress(ss_sector.str().c_str(), &a_defaultSectorY[iSector]);
          defaultTreeY->GetEntry(0);
        }
        sectorNameTree->SetBranchAddress(ss_sector.str().c_str(), &a_sectorName[iSector]);
        sectorNameTree->GetEntry(0);
      }
    }

    // Check whether sector definitions are identical with the ones of previous iterations and with the ones in baseline file
    for (auto& i_sector : m_tkSector_) {
      const std::string& name(i_sector.second.name);
      if (!firstIter) {
        const std::string& nameLastIter(*a_sectorName[i_sector.first]);
        if (name != nameLastIter) {
          edm::LogError("CalculateAPE") << "Inconsistent sector definition in iterationFile for sector "
                                        << i_sector.first << ",\n"
                                        << "Recent iteration has name \"" << name << "\", while previous had \""
                                        << nameLastIter << "\"\n"
                                        << "...APE calculation stopped. Please check sector definitions in config!\n";
          return;
        }
      } else
        a_sectorName[i_sector.first] = new std::string(name);
      if (!setBaseline && baselineFile) {
        const std::string& nameBaseline(*a_sectorBaselineName[i_sector.first]);
        if (name != nameBaseline) {
          edm::LogError("CalculateAPE") << "Inconsistent sector definition in baselineFile for sector "
                                        << i_sector.first << ",\n"
                                        << "Recent iteration has name \"" << name << "\", while baseline had \""
                                        << nameBaseline << "\"\n"
                                        << "...APE calculation stopped. Please check sector definitions in config!\n";
          return;
        }
      }
    }

    // Loop over sectors for calculating getting default APE
    for (auto& i_sector : m_tkSector_) {
      double defaultApeX(0.);
      double defaultApeY(0.);
      unsigned int nModules(0);
      for (auto const& i_rawId : i_sector.second.v_rawId) {
        std::vector<AlignTransformErrorExtended> alignErrors = alignmentErrors->m_alignError;
        for (std::vector<AlignTransformErrorExtended>::const_iterator i_alignError = alignErrors.begin();
             i_alignError != alignErrors.end();
             ++i_alignError) {
          if (i_rawId == i_alignError->rawId()) {
            CLHEP::HepSymMatrix errMatrix = i_alignError->matrix();
            defaultApeX += errMatrix[0][0];
            defaultApeY += errMatrix[1][1];
            nModules++;
          }
        }
      }
      a_defaultSectorX[i_sector.first] = defaultApeX / nModules;
      a_defaultSectorY[i_sector.first] = defaultApeY / nModules;
    }
    sectorNameTree->Fill();
    sectorNameTree->Write("nameTree");
    defaultTreeX->Fill();
    defaultTreeX->Write("iterTreeX");
    defaultTreeY->Fill();
    defaultTreeY->Write("iterTreeY");

    defaultFile->Close();
    delete defaultFile;

    if (baselineFile) {
      baselineFile->Close();
      delete baselineFile;
    }

    firstEvent = false;
  }
}

// ------------ method called once each job just before starting event loop  ------------
void ApeEstimatorSummary::beginJob() {
  this->openInputFile();

  this->getTrackerSectorStructs();

  this->bookHists();
}

// ------------ method called once each job just after ending the event loop  ------------
void ApeEstimatorSummary::endJob() {
  this->calculateApe();

  this->writeHists();
}

//define this as a plug-in
DEFINE_FWK_MODULE(ApeEstimatorSummary);