Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 // -*- C++ -*-
0002 //
0003 // Package:    ApeEstimator
0004 // Class:      ApeEstimator
0005 //
0006 /**\class ApeEstimator ApeEstimator.cc Alignment/APEEstimation/src/ApeEstimator.cc
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Johannes Hauk
0015 //         Created:  Tue Jan  6 15:02:09 CET 2009
0016 //         Modified by: Christian Schomakers (RWTH Aachen)
0017 // $Id: ApeEstimator.cc,v 1.27 2012/06/26 09:42:33 hauk Exp $
0018 //
0019 //
0020 
0021 // system include files
0022 #include <memory>
0023 #include <sstream>
0024 #include <fstream>
0025 
0026 // user include files
0027 #include "FWCore/Framework/interface/Frameworkfwd.h"
0028 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/EventSetup.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0034 #include "FWCore/ServiceRegistry/interface/Service.h"
0035 #include "FWCore/Utilities/interface/InputTag.h"
0036 #include "FWCore/Utilities/interface/EDGetToken.h"
0037 
0038 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0039 #include "CommonTools/Utils/interface/TFileDirectory.h"
0040 
0041 #include "DataFormats/Common/interface/Handle.h"
0042 #include "DataFormats/TrackReco/interface/Track.h"
0043 #include "DataFormats/TrackReco/interface/HitPattern.h"
0044 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0045 #include "DataFormats/DetId/interface/DetId.h"
0046 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0047 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0048 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0049 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0050 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0051 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0052 #include "DataFormats/GeometrySurface/interface/LocalError.h"
0053 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
0054 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h"
0055 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0056 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0057 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0058 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0059 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0060 
0061 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0062 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0063 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0064 #include "TrackingTools/PatternTools/interface/MeasurementExtractor.h"
0065 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0066 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0067 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0068 
0069 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0070 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0071 #include "Geometry/CommonTopologies/interface/RadialStripTopology.h"
0072 //added by Ajay 6Nov 2014
0073 //.......................
0074 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0075 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0076 
0077 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorExtendedRcd.h"
0078 
0079 //...............
0080 //
0081 
0082 ////.....................
0083 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0084 
0085 #include "CondFormats/Alignment/interface/Definitions.h"
0086 
0087 #include "RecoLocalTracker/SiStripClusterizer/interface/SiStripClusterInfo.h"
0088 
0089 #include "Alignment/APEEstimation/interface/TrackerSectorStruct.h"
0090 #include "Alignment/APEEstimation/interface/TrackerDetectorStruct.h"
0091 #include "Alignment/APEEstimation/interface/EventVariables.h"
0092 #include "Alignment/APEEstimation/interface/ReducedTrackerTreeVariables.h"
0093 
0094 #include "TH1.h"
0095 #include "TH2.h"
0096 #include "TProfile.h"
0097 #include "TFile.h"
0098 #include "TTree.h"
0099 #include "TF1.h"
0100 #include "TString.h"
0101 #include "TMath.h"
0102 
0103 #include "RecoLocalTracker/SiStripRecHitConverter/interface/StripCPE.h"
0104 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0105 
0106 #include "MagneticField/Engine/interface/MagneticField.h"
0107 #include "CondFormats/SiStripObjects/interface/SiStripLorentzAngle.h"
0108 //ADDED BY LOIC QUERTENMONT
0109 #include "FWCore/Framework/interface/DependentRecordImplementation.h"
0110 #include "CalibTracker/Records/interface/SiStripDependentRecords.h"
0111 /////////
0112 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0113 #include "DataFormats/GeometrySurface/interface/Bounds.h"
0114 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0115 #include "CondFormats/DataRecord/interface/SiStripLorentzAngleRcd.h"
0116 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0117 #include "DataFormats/Math/interface/Point3D.h"
0118 #include "FWCore/Utilities/interface/EDMException.h"
0119 
0120 //
0121 // class decleration
0122 //
0123 
0124 class ApeEstimator : public edm::one::EDAnalyzer<> {
0125 public:
0126   explicit ApeEstimator(const edm::ParameterSet&);
0127   ~ApeEstimator() override;
0128 
0129 private:
0130   struct PositionAndError2 {
0131     PositionAndError2() : posX(-999.F), posY(-999.F), errX2(-999.F), errY2(-999.F){};
0132     PositionAndError2(float x, float y, float eX, float eY) : posX(x), posY(y), errX2(eX), errY2(eY){};
0133     float posX;
0134     float posY;
0135     float errX2;
0136     float errY2;
0137   };
0138   typedef std::pair<TrackStruct::HitState, PositionAndError2> StatePositionAndError2;
0139 
0140   void beginJob() override;
0141   void analyze(const edm::Event&, const edm::EventSetup&) override;
0142   void endJob() override;
0143 
0144   bool isHit2D(const TrackingRecHit&) const;
0145 
0146   void sectorBuilder();
0147   bool checkIntervalsForSectors(const unsigned int sectorCounter, const std::vector<double>&) const;
0148   bool checkModuleIds(const unsigned int, const std::vector<unsigned int>&) const;
0149   bool checkModuleBools(const bool, const std::vector<unsigned int>&) const;
0150   bool checkModuleDirections(const int, const std::vector<int>&) const;
0151   bool checkModulePositions(const float, const std::vector<double>&) const;
0152   void statistics(const TrackerSectorStruct&, const int) const;
0153 
0154   void residualErrorBinning();
0155 
0156   void bookSectorHistsForAnalyzerMode();
0157   void bookSectorHistsForApeCalculation();
0158   void bookTrackHists();
0159 
0160   TrackStruct::TrackParameterStruct fillTrackVariables(const reco::Track&, const Trajectory&, const reco::BeamSpot&);
0161   TrackStruct::HitParameterStruct fillHitVariables(const TrajectoryMeasurement&, const edm::EventSetup&);
0162 
0163   StatePositionAndError2 positionAndError2(const LocalPoint&, const LocalError&, const TransientTrackingRecHit&);
0164   PositionAndError2 rectangularPositionAndError2(const LocalPoint&, const LocalError&);
0165   PositionAndError2 radialPositionAndError2(const LocalPoint&, const LocalError&, const RadialStripTopology&);
0166 
0167   void hitSelection();
0168   void setHitSelectionMap(const std::string&);
0169   void setHitSelectionMapUInt(const std::string&);
0170   bool hitSelected(TrackStruct::HitParameterStruct&) const;
0171   bool inDoubleInterval(const std::vector<double>&, const float) const;
0172   bool inUintInterval(const std::vector<unsigned int>&, const unsigned int, const unsigned int = 999) const;
0173 
0174   void fillHistsForAnalyzerMode(const TrackStruct&);
0175   void fillHitHistsXForAnalyzerMode(const TrackStruct::HitParameterStruct&, TrackerSectorStruct&);
0176   void fillHitHistsYForAnalyzerMode(const TrackStruct::HitParameterStruct&, TrackerSectorStruct&);
0177   void fillHistsForApeCalculation(const TrackStruct&);
0178 
0179   void calculateAPE();
0180 
0181   // ----------member data ---------------------------
0182   const edm::ParameterSet parameterSet_;
0183   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0184   const edm::ESGetToken<SiStripLorentzAngle, SiStripLorentzAngleDepRcd> lorentzAngleToken_;
0185 
0186   std::map<unsigned int, TrackerSectorStruct> m_tkSector_;
0187   TrackerDetectorStruct tkDetector_;
0188   SiStripClusterInfo siStripClusterInfo_;
0189 
0190   edm::EDGetTokenT<TrajTrackAssociationCollection> tjTagToken_;
0191   edm::EDGetTokenT<reco::BeamSpot> offlinebeamSpot_;
0192 
0193   std::map<unsigned int, std::pair<double, double> > m_resErrBins_;
0194   std::map<unsigned int, ReducedTrackerTreeVariables> m_tkTreeVar_;
0195 
0196   std::map<std::string, std::vector<double> > m_hitSelection_;
0197   std::map<std::string, std::vector<unsigned int> > m_hitSelectionUInt_;
0198 
0199   bool trackCut_;
0200 
0201   const unsigned int maxTracksPerEvent_;
0202   const unsigned int minGoodHitsPerTrack_;
0203 
0204   const bool analyzerMode_;
0205 
0206   const bool calculateApe_;
0207 
0208   unsigned int counter1, counter2, counter3, counter4, counter5, counter6;
0209 };
0210 
0211 //
0212 // constants, enums and typedefs
0213 //
0214 
0215 //
0216 // static data member definitions
0217 //
0218 
0219 //
0220 // constructors and destructor
0221 //
0222 ApeEstimator::ApeEstimator(const edm::ParameterSet& iConfig)
0223     : parameterSet_(iConfig),
0224       magFieldToken_(esConsumes()),
0225       lorentzAngleToken_(esConsumes()),
0226       siStripClusterInfo_(consumesCollector(), std::string("")),
0227       tjTagToken_(
0228           consumes<TrajTrackAssociationCollection>(parameterSet_.getParameter<edm::InputTag>("tjTkAssociationMapTag"))),
0229       offlinebeamSpot_(consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"))),
0230       trackCut_(false),
0231       maxTracksPerEvent_(parameterSet_.getParameter<unsigned int>("maxTracksPerEvent")),
0232       minGoodHitsPerTrack_(parameterSet_.getParameter<unsigned int>("minGoodHitsPerTrack")),
0233       analyzerMode_(parameterSet_.getParameter<bool>("analyzerMode")),
0234       calculateApe_(parameterSet_.getParameter<bool>("calculateApe")) {
0235   counter1 = counter2 = counter3 = counter4 = counter5 = counter6 = 0;
0236 }
0237 
0238 ApeEstimator::~ApeEstimator() {}
0239 
0240 //
0241 // member functions
0242 //
0243 
0244 // -----------------------------------------------------------------------------------------------------------
0245 
0246 void ApeEstimator::sectorBuilder() {
0247   TFile* tkTreeFile(TFile::Open((parameterSet_.getParameter<std::string>("TrackerTreeFile")).c_str()));
0248   if (tkTreeFile) {
0249     edm::LogInfo("SectorBuilder") << "TrackerTreeFile OK";
0250   } else {
0251     edm::LogError("SectorBuilder") << "TrackerTreeFile not found";
0252     return;
0253   }
0254   TTree* tkTree(nullptr);
0255   tkTreeFile->GetObject("TrackerTreeGenerator/TrackerTree/TrackerTree", tkTree);
0256   if (tkTree) {
0257     edm::LogInfo("SectorBuilder") << "TrackerTree OK";
0258   } else {
0259     edm::LogError("SectorBuilder") << "TrackerTree not found in file";
0260     return;
0261   }
0262   unsigned int rawId(999), subdetId(999), layer(999), side(999), half(999), rod(999), ring(999), petal(999), blade(999),
0263       panel(999), outerInner(999), module(999), nStrips(999);
0264   bool isDoubleSide(false), isRPhi(false), isStereo(false);
0265   int uDirection(999), vDirection(999), wDirection(999);
0266   float posR(999.F), posPhi(999.F), posEta(999.F), posX(999.F), posY(999.F), posZ(999.F);
0267 
0268   tkTree->SetBranchAddress("RawId", &rawId);
0269   tkTree->SetBranchAddress("SubdetId", &subdetId);
0270   tkTree->SetBranchAddress("Layer", &layer);
0271   tkTree->SetBranchAddress("Side", &side);
0272   tkTree->SetBranchAddress("Half", &half);
0273   tkTree->SetBranchAddress("Rod", &rod);
0274   tkTree->SetBranchAddress("Ring", &ring);
0275   tkTree->SetBranchAddress("Petal", &petal);
0276   tkTree->SetBranchAddress("Blade", &blade);
0277   tkTree->SetBranchAddress("Panel", &panel);
0278   tkTree->SetBranchAddress("OuterInner", &outerInner);
0279   tkTree->SetBranchAddress("Module", &module);
0280   tkTree->SetBranchAddress("NStrips", &nStrips);
0281   tkTree->SetBranchAddress("IsDoubleSide", &isDoubleSide);
0282   tkTree->SetBranchAddress("IsRPhi", &isRPhi);
0283   tkTree->SetBranchAddress("IsStereo", &isStereo);
0284   tkTree->SetBranchAddress("UDirection", &uDirection);
0285   tkTree->SetBranchAddress("VDirection", &vDirection);
0286   tkTree->SetBranchAddress("WDirection", &wDirection);
0287   tkTree->SetBranchAddress("PosR", &posR);
0288   tkTree->SetBranchAddress("PosPhi", &posPhi);
0289   tkTree->SetBranchAddress("PosEta", &posEta);
0290   tkTree->SetBranchAddress("PosX", &posX);
0291   tkTree->SetBranchAddress("PosY", &posY);
0292   tkTree->SetBranchAddress("PosZ", &posZ);
0293 
0294   int nModules(tkTree->GetEntries());
0295   TrackerSectorStruct allSectors;
0296 
0297   //Loop over all Sectors
0298   unsigned int sectorCounter(0);
0299   std::vector<edm::ParameterSet> v_sectorDef(parameterSet_.getParameter<std::vector<edm::ParameterSet> >("Sectors"));
0300   edm::LogInfo("SectorBuilder") << "There are " << v_sectorDef.size() << " Sectors definded";
0301   for (auto const& parSet : v_sectorDef) {
0302     ++sectorCounter;
0303     const std::string& sectorName(parSet.getParameter<std::string>("name"));
0304     std::vector<unsigned int> v_rawId(parSet.getParameter<std::vector<unsigned int> >("rawId")),
0305         v_subdetId(parSet.getParameter<std::vector<unsigned int> >("subdetId")),
0306         v_layer(parSet.getParameter<std::vector<unsigned int> >("layer")),
0307         v_side(parSet.getParameter<std::vector<unsigned int> >("side")),
0308         v_half(parSet.getParameter<std::vector<unsigned int> >("half")),
0309         v_rod(parSet.getParameter<std::vector<unsigned int> >("rod")),
0310         v_ring(parSet.getParameter<std::vector<unsigned int> >("ring")),
0311         v_petal(parSet.getParameter<std::vector<unsigned int> >("petal")),
0312         v_blade(parSet.getParameter<std::vector<unsigned int> >("blade")),
0313         v_panel(parSet.getParameter<std::vector<unsigned int> >("panel")),
0314         v_outerInner(parSet.getParameter<std::vector<unsigned int> >("outerInner")),
0315         v_module(parSet.getParameter<std::vector<unsigned int> >("module")),
0316         v_nStrips(parSet.getParameter<std::vector<unsigned int> >("nStrips")),
0317         v_isDoubleSide(parSet.getParameter<std::vector<unsigned int> >("isDoubleSide")),
0318         v_isRPhi(parSet.getParameter<std::vector<unsigned int> >("isRPhi")),
0319         v_isStereo(parSet.getParameter<std::vector<unsigned int> >("isStereo"));
0320     std::vector<int> v_uDirection(parSet.getParameter<std::vector<int> >("uDirection")),
0321         v_vDirection(parSet.getParameter<std::vector<int> >("vDirection")),
0322         v_wDirection(parSet.getParameter<std::vector<int> >("wDirection"));
0323     std::vector<double> v_posR(parSet.getParameter<std::vector<double> >("posR")),
0324         v_posPhi(parSet.getParameter<std::vector<double> >("posPhi")),
0325         v_posEta(parSet.getParameter<std::vector<double> >("posEta")),
0326         v_posX(parSet.getParameter<std::vector<double> >("posX")),
0327         v_posY(parSet.getParameter<std::vector<double> >("posY")),
0328         v_posZ(parSet.getParameter<std::vector<double> >("posZ"));
0329 
0330     if (!this->checkIntervalsForSectors(sectorCounter, v_posR) ||
0331         !this->checkIntervalsForSectors(sectorCounter, v_posPhi) ||
0332         !this->checkIntervalsForSectors(sectorCounter, v_posEta) ||
0333         !this->checkIntervalsForSectors(sectorCounter, v_posX) ||
0334         !this->checkIntervalsForSectors(sectorCounter, v_posY) ||
0335         !this->checkIntervalsForSectors(sectorCounter, v_posZ))
0336       continue;
0337 
0338     TrackerSectorStruct tkSector;
0339     tkSector.name = sectorName;
0340 
0341     ReducedTrackerTreeVariables tkTreeVar;
0342 
0343     //Loop over all Modules
0344     for (int module = 0; module < nModules; ++module) {
0345       tkTree->GetEntry(module);
0346 
0347       if (sectorCounter == 1) {
0348         tkTreeVar.subdetId = subdetId;
0349         tkTreeVar.nStrips = nStrips;
0350         tkTreeVar.uDirection = uDirection;
0351         tkTreeVar.vDirection = vDirection;
0352         tkTreeVar.wDirection = wDirection;
0353         m_tkTreeVar_[rawId] = tkTreeVar;
0354       }
0355 
0356       if (!this->checkModuleIds(rawId, v_rawId))
0357         continue;
0358       if (!this->checkModuleIds(subdetId, v_subdetId))
0359         continue;
0360       if (!this->checkModuleIds(layer, v_layer))
0361         continue;
0362       if (!this->checkModuleIds(side, v_side))
0363         continue;
0364       if (!this->checkModuleIds(half, v_half))
0365         continue;
0366       if (!this->checkModuleIds(rod, v_rod))
0367         continue;
0368       if (!this->checkModuleIds(ring, v_ring))
0369         continue;
0370       if (!this->checkModuleIds(petal, v_petal))
0371         continue;
0372       if (!this->checkModuleIds(blade, v_blade))
0373         continue;
0374       if (!this->checkModuleIds(panel, v_panel))
0375         continue;
0376       if (!this->checkModuleIds(outerInner, v_outerInner))
0377         continue;
0378       if (!this->checkModuleIds(module, v_module))
0379         continue;
0380       if (!this->checkModuleIds(nStrips, v_nStrips))
0381         continue;
0382       if (!this->checkModuleBools(isDoubleSide, v_isDoubleSide))
0383         continue;
0384       if (!this->checkModuleBools(isRPhi, v_isRPhi))
0385         continue;
0386       if (!this->checkModuleBools(isStereo, v_isStereo))
0387         continue;
0388       if (!this->checkModuleDirections(uDirection, v_uDirection))
0389         continue;
0390       if (!this->checkModuleDirections(vDirection, v_vDirection))
0391         continue;
0392       if (!this->checkModuleDirections(wDirection, v_wDirection))
0393         continue;
0394       if (!this->checkModulePositions(posR, v_posR))
0395         continue;
0396       if (!this->checkModulePositions(posPhi, v_posPhi))
0397         continue;
0398       if (!this->checkModulePositions(posEta, v_posEta))
0399         continue;
0400       if (!this->checkModulePositions(posX, v_posX))
0401         continue;
0402       if (!this->checkModulePositions(posY, v_posY))
0403         continue;
0404       if (!this->checkModulePositions(posZ, v_posZ))
0405         continue;
0406 
0407       tkSector.v_rawId.push_back(rawId);
0408       bool moduleSelected(false);
0409       for (auto const& i_rawId : allSectors.v_rawId) {
0410         if (rawId == i_rawId)
0411           moduleSelected = true;
0412       }
0413       if (!moduleSelected)
0414         allSectors.v_rawId.push_back(rawId);
0415     }
0416 
0417     bool isPixel(false);
0418     bool isStrip(false);
0419     for (auto const& i_rawId : tkSector.v_rawId) {
0420       switch (m_tkTreeVar_[i_rawId].subdetId) {
0421         case PixelSubdetector::PixelBarrel:
0422         case PixelSubdetector::PixelEndcap:
0423           isPixel = true;
0424           break;
0425         case StripSubdetector::TIB:
0426         case StripSubdetector::TOB:
0427         case StripSubdetector::TID:
0428         case StripSubdetector::TEC:
0429           isStrip = true;
0430           break;
0431       }
0432     }
0433 
0434     if (isPixel && isStrip) {
0435       edm::LogError("SectorBuilder")
0436           << "Incorrect Sector Definition: there are pixel and strip modules within one sector"
0437           << "\n... sector selection is not applied, sector " << sectorCounter << " is not built";
0438       continue;
0439     }
0440     tkSector.isPixel = isPixel;
0441 
0442     m_tkSector_[sectorCounter] = tkSector;
0443     edm::LogInfo("SectorBuilder") << "There are " << tkSector.v_rawId.size() << " Modules in Sector " << sectorCounter;
0444   }
0445   this->statistics(allSectors, nModules);
0446   return;
0447 }
0448 
0449 // -----------------------------------------------------------------------------------------------------------
0450 
0451 bool ApeEstimator::checkIntervalsForSectors(const unsigned int sectorCounter, const std::vector<double>& v_id) const {
0452   if (v_id.empty())
0453     return true;
0454   if (v_id.size() % 2 == 1) {
0455     edm::LogError("SectorBuilder")
0456         << "Incorrect Sector Definition: Position Vectors need even number of arguments (Intervals)"
0457         << "\n... sector selection is not applied, sector " << sectorCounter << " is not built";
0458     return false;
0459   }
0460   int entry(0);
0461   double intervalBegin(999.);
0462   for (auto const& i_id : v_id) {
0463     ++entry;
0464     if (entry % 2 == 1)
0465       intervalBegin = i_id;
0466     if (entry % 2 == 0 && intervalBegin > i_id) {
0467       edm::LogError("SectorBuilder") << "Incorrect Sector Definition (Position Vector Intervals): \t" << intervalBegin
0468                                      << " is bigger than " << i_id << " but is expected to be smaller"
0469                                      << "\n... sector selection is not applied, sector " << sectorCounter
0470                                      << " is not built";
0471       return false;
0472     }
0473   }
0474   return true;
0475 }
0476 
0477 bool ApeEstimator::checkModuleIds(const unsigned int id, const std::vector<unsigned int>& v_id) const {
0478   if (v_id.empty())
0479     return true;
0480   for (auto const& i_id : v_id) {
0481     if (id == i_id)
0482       return true;
0483   }
0484   return false;
0485 }
0486 
0487 bool ApeEstimator::checkModuleBools(const bool id, const std::vector<unsigned int>& v_id) const {
0488   if (v_id.empty())
0489     return true;
0490   for (auto const& i_id : v_id) {
0491     if (1 == i_id && id)
0492       return true;
0493     if (2 == i_id && !id)
0494       return true;
0495   }
0496   return false;
0497 }
0498 
0499 bool ApeEstimator::checkModuleDirections(const int id, const std::vector<int>& v_id) const {
0500   if (v_id.empty())
0501     return true;
0502   for (auto const& i_id : v_id) {
0503     if (id == i_id)
0504       return true;
0505   }
0506   return false;
0507 }
0508 
0509 bool ApeEstimator::checkModulePositions(const float id, const std::vector<double>& v_id) const {
0510   if (v_id.empty())
0511     return true;
0512   int entry(0);
0513   double intervalBegin(999.);
0514   for (auto const& i_id : v_id) {
0515     ++entry;
0516     if (entry % 2 == 1)
0517       intervalBegin = i_id;
0518     if (entry % 2 == 0 && id >= intervalBegin && id < i_id)
0519       return true;
0520   }
0521   return false;
0522 }
0523 
0524 void ApeEstimator::statistics(const TrackerSectorStruct& allSectors, const int nModules) const {
0525   bool commonModules(false);
0526   for (std::map<unsigned int, TrackerSectorStruct>::const_iterator i_sector = m_tkSector_.begin();
0527        i_sector != m_tkSector_.end();
0528        ++i_sector) {
0529     std::map<unsigned int, TrackerSectorStruct>::const_iterator i_sector2(i_sector);
0530     for (++i_sector2; i_sector2 != m_tkSector_.end(); ++i_sector2) {
0531       unsigned int nCommonModules(0);
0532       for (auto const& i_module : (*i_sector).second.v_rawId) {
0533         for (auto const& i_module2 : (*i_sector2).second.v_rawId) {
0534           if (i_module2 == i_module)
0535             ++nCommonModules;
0536         }
0537       }
0538       if (nCommonModules == 0)
0539         ;  //edm::LogInfo("SectorBuilder")<<"Sector "<<(*i_sector).first<<" and Sector "<<(*i_sector2).first<< " have ZERO Modules in common";
0540       else {
0541         edm::LogError("SectorBuilder") << "Sector " << (*i_sector).first << " and Sector " << (*i_sector2).first
0542                                        << " have " << nCommonModules << " Modules in common";
0543         commonModules = true;
0544       }
0545     }
0546   }
0547   if (static_cast<int>(allSectors.v_rawId.size()) == nModules)
0548     edm::LogInfo("SectorBuilder") << "ALL Tracker Modules are contained in the Sectors";
0549   else
0550     edm::LogWarning("SectorBuilder") << "There are " << allSectors.v_rawId.size() << " Modules in all Sectors"
0551                                      << " out of " << nModules << " Tracker Modules";
0552   if (!commonModules)
0553     edm::LogInfo("SectorBuilder") << "There are ZERO modules associated to different sectors, no ambiguities exist";
0554   else
0555     edm::LogError("SectorBuilder")
0556         << "There are modules associated to different sectors, APE value cannot be assigned reasonably";
0557 }
0558 
0559 // -----------------------------------------------------------------------------------------------------------
0560 
0561 void ApeEstimator::residualErrorBinning() {
0562   std::vector<double> v_residualErrorBinning(parameterSet_.getParameter<std::vector<double> >("residualErrorBinning"));
0563   if (v_residualErrorBinning.size() == 1) {
0564     edm::LogError("ResidualErrorBinning") << "Incorrect selection of Residual Error Bins (used for APE calculation): \t"
0565                                           << "Only one argument passed, so no interval is specified"
0566                                           << "\n... delete whole bin selection";  //m_resErrBins_ remains empty
0567     return;
0568   }
0569   double xMin(0.), xMax(0.);
0570   unsigned int binCounter(-1);
0571   for (auto const& i_binning : v_residualErrorBinning) {
0572     ++binCounter;
0573     if (binCounter == 0) {
0574       xMin = i_binning;
0575       continue;
0576     }
0577     xMax = i_binning;
0578     if (xMax <= xMin) {
0579       edm::LogError("ResidualErrorBinning")
0580           << "Incorrect selection of Residual Error Bins (used for APE calculation): \t" << xMin << " is bigger than "
0581           << xMax << " but is expected to be smaller"
0582           << "\n... delete whole bin selection";
0583       m_resErrBins_.clear();
0584       return;
0585     }
0586     m_resErrBins_[binCounter].first = xMin;
0587     m_resErrBins_[binCounter].second = xMax;
0588     xMin = xMax;
0589   }
0590   edm::LogInfo("ResidualErrorBinning")
0591       << m_resErrBins_.size() << " Intervals of residual errors used for separate APE calculation sucessfully set";
0592 }
0593 
0594 // -----------------------------------------------------------------------------------------------------------
0595 
0596 void ApeEstimator::bookSectorHistsForAnalyzerMode() {
0597   std::vector<unsigned int> v_errHists(parameterSet_.getParameter<std::vector<unsigned int> >("vErrHists"));
0598   for (std::vector<unsigned int>::iterator i_errHists = v_errHists.begin(); i_errHists != v_errHists.end();
0599        ++i_errHists) {
0600     for (std::vector<unsigned int>::iterator i_errHists2 = i_errHists; i_errHists2 != v_errHists.end();) {
0601       ++i_errHists2;
0602       if (*i_errHists == *i_errHists2) {
0603         edm::LogError("BookSectorHists") << "Value of vErrHists in config exists twice: " << *i_errHists
0604                                          << "\n... delete one of both";
0605         v_errHists.erase(i_errHists2);
0606       }
0607     }
0608   }
0609 
0610   for (auto& i_sector : m_tkSector_) {
0611     bool zoomHists(parameterSet_.getParameter<bool>("zoomHists"));
0612 
0613     double widthMax = zoomHists ? 20. : 200.;
0614     double chargePixelMax = zoomHists ? 200000. : 2000000.;
0615     double chargeStripMax = zoomHists ? 1000. : 10000.;
0616     double sOverNMax = zoomHists ? 200. : 2000.;
0617     double logClusterProbMin = zoomHists ? -5. : -15.;
0618 
0619     double resXAbsMax = zoomHists ? 0.5 : 5.;
0620     double norResXAbsMax = zoomHists ? 10. : 50.;
0621     double probXMin = zoomHists ? -0.01 : -0.1;
0622     double probXMax = zoomHists ? 0.11 : 1.1;
0623     double sigmaXMin = zoomHists ? 0. : -0.05;
0624     double sigmaXMax = zoomHists ? 0.02 : 1.;
0625     double sigmaX2Max = sigmaXMax * sigmaXMax;
0626     double sigmaXHitMax = zoomHists ? 0.02 : 1.;
0627     double phiSensXMax = zoomHists ? 31. : 93.;
0628 
0629     double norChi2Max = zoomHists ? 5. : 1000.;
0630     double d0Max = zoomHists ? 0.02 : 40.;  // cosmics: 100.|100.
0631     double dzMax = zoomHists ? 15. : 100.;  // cosmics: 200.|600.
0632     double pMax = zoomHists ? 200. : 2000.;
0633     double invPMax = zoomHists ? 0.05 : 10.;  //begins at 20GeV, 0.1GeV
0634 
0635     edm::Service<TFileService> fileService;
0636     if (!fileService) {
0637       throw edm::Exception(edm::errors::Configuration, "TFileService is not registered in cfg file");
0638     }
0639 
0640     std::stringstream sector;
0641     sector << "Sector_" << i_sector.first;
0642     TFileDirectory secDir = fileService->mkdir(sector.str());
0643 
0644     // Dummy histo containing the sector name as title
0645     i_sector.second.Name = secDir.make<TH1F>("z_name", i_sector.second.name.c_str(), 1, 0, 1);
0646 
0647     // Do not book histos for empty sectors
0648     if (i_sector.second.v_rawId.empty()) {
0649       continue;
0650     }
0651     // Set parameters for correlationHists
0652     i_sector.second.setCorrHistParams(&secDir, norResXAbsMax, sigmaXHitMax, sigmaXMax);
0653 
0654     // Book pixel or strip specific hists
0655     const bool pixelSector(i_sector.second.isPixel);
0656 
0657     // Cluster Parameters
0658     i_sector.second.m_correlationHistsX["WidthX"] = i_sector.second.bookCorrHistsX("WidthX",
0659                                                                                    "cluster width",
0660                                                                                    "w_{cl,x}",
0661                                                                                    "[# channels]",
0662                                                                                    static_cast<int>(widthMax),
0663                                                                                    static_cast<int>(widthMax),
0664                                                                                    0.,
0665                                                                                    widthMax,
0666                                                                                    "nph");
0667     i_sector.second.m_correlationHistsX["BaryStripX"] = i_sector.second.bookCorrHistsX(
0668         "BaryStripX", "barycenter of cluster charge", "b_{cl,x}", "[# channels]", 800, 100, -10., 790., "nph");
0669 
0670     if (pixelSector) {
0671       i_sector.second.m_correlationHistsY["WidthY"] = i_sector.second.bookCorrHistsY("WidthY",
0672                                                                                      "cluster width",
0673                                                                                      "w_{cl,y}",
0674                                                                                      "[# channels]",
0675                                                                                      static_cast<int>(widthMax),
0676                                                                                      static_cast<int>(widthMax),
0677                                                                                      0.,
0678                                                                                      widthMax,
0679                                                                                      "nph");
0680       i_sector.second.m_correlationHistsY["BaryStripY"] = i_sector.second.bookCorrHistsY(
0681           "BaryStripY", "barycenter of cluster charge", "b_{cl,y}", "[# channels]", 800, 100, -10., 790., "nph");
0682 
0683       i_sector.second.m_correlationHistsX["ChargePixel"] = i_sector.second.bookCorrHistsX(
0684           "ChargePixel", "cluster charge", "c_{cl}", "[e]", 100, 50, 0., chargePixelMax, "nph");
0685       i_sector.second.m_correlationHistsX["ClusterProbXY"] = i_sector.second.bookCorrHistsX(
0686           "ClusterProbXY", "cluster probability xy", "prob_{cl,xy}", "", 100, 50, 0., 1., "nph");
0687       i_sector.second.m_correlationHistsX["ClusterProbQ"] = i_sector.second.bookCorrHistsX(
0688           "ClusterProbQ", "cluster probability q", "prob_{cl,q}", "", 100, 50, 0., 1., "nph");
0689       i_sector.second.m_correlationHistsX["ClusterProbXYQ"] = i_sector.second.bookCorrHistsX(
0690           "ClusterProbXYQ", "cluster probability xyq", "prob_{cl,xyq}", "", 100, 50, 0., 1., "nph");
0691       i_sector.second.m_correlationHistsX["LogClusterProb"] = i_sector.second.bookCorrHistsX(
0692           "LogClusterProb", "cluster probability xy", "log(prob_{cl,xy})", "", 60, 30, logClusterProbMin, 0., "nph");
0693       i_sector.second.m_correlationHistsX["IsOnEdge"] =
0694           i_sector.second.bookCorrHistsX("IsOnEdge", "IsOnEdge", "isOnEdge", "", 2, 2, 0, 2, "nph");
0695       i_sector.second.m_correlationHistsX["HasBadPixels"] =
0696           i_sector.second.bookCorrHistsX("HasBadPixels", "HasBadPixels", "hasBadPixels", "", 2, 2, 0, 2, "nph");
0697       i_sector.second.m_correlationHistsX["SpansTwoRoc"] =
0698           i_sector.second.bookCorrHistsX("SpansTwoRoc", "SpansTwoRoc", "spansTwoRoc", "", 2, 2, 0, 2, "nph");
0699       i_sector.second.m_correlationHistsX["QBin"] =
0700           i_sector.second.bookCorrHistsX("QBin", "q bin", "q bin", "", 8, 8, 0, 8, "nph");
0701 
0702       i_sector.second.m_correlationHistsY["ChargePixel"] = i_sector.second.bookCorrHistsY(
0703           "ChargePixel", "cluster charge", "c_{cl}", "[e]", 100, 50, 0., chargePixelMax, "nph");
0704       i_sector.second.m_correlationHistsY["ClusterProbXY"] = i_sector.second.bookCorrHistsY(
0705           "ClusterProbXY", "cluster probability xy", "prob_{cl,xy}", "", 100, 50, 0., 1., "nph");
0706       i_sector.second.m_correlationHistsY["ClusterProbQ"] = i_sector.second.bookCorrHistsY(
0707           "ClusterProbQ", "cluster probability q", "prob_{cl,q}", "", 100, 50, 0., 1., "nph");
0708       i_sector.second.m_correlationHistsY["ClusterProbXYQ"] = i_sector.second.bookCorrHistsY(
0709           "ClusterProbXYQ", "cluster probability xyq", "prob_{cl,xyq}", "", 100, 50, 0., 1., "nph");
0710       i_sector.second.m_correlationHistsY["LogClusterProb"] = i_sector.second.bookCorrHistsY(
0711           "LogClusterProb", "cluster probability xy", "log(prob_{cl,xy})", "", 60, 30, logClusterProbMin, 0., "nph");
0712       i_sector.second.m_correlationHistsY["IsOnEdge"] =
0713           i_sector.second.bookCorrHistsY("IsOnEdge", "IsOnEdge", "isOnEdge", "", 2, 2, 0, 2, "nph");
0714       i_sector.second.m_correlationHistsY["HasBadPixels"] =
0715           i_sector.second.bookCorrHistsY("HasBadPixels", "HasBadPixels", "hasBadPixels", "", 2, 2, 0, 2, "nph");
0716       i_sector.second.m_correlationHistsY["SpansTwoRoc"] =
0717           i_sector.second.bookCorrHistsY("SpansTwoRoc", "SpansTwoRoc", "spansTwoRoc", "", 2, 2, 0, 2, "nph");
0718       i_sector.second.m_correlationHistsY["QBin"] =
0719           i_sector.second.bookCorrHistsY("QBin", "q bin", "q bin", "", 8, 8, 0, 8, "nph");
0720     }
0721 
0722     else {
0723       i_sector.second.m_correlationHistsX["ChargeStrip"] = i_sector.second.bookCorrHistsX(
0724           "ChargeStrip", "cluster charge", "c_{cl}", "[APV counts]", 100, 50, 0., chargeStripMax, "nph");
0725       i_sector.second.m_correlationHistsX["MaxStrip"] = i_sector.second.bookCorrHistsX(
0726           "MaxStrip", "strip with max. charge", "n_{cl,max}", "[# strips]", 800, 800, -10., 790., "npht");
0727       i_sector.second.m_correlationHistsX["MaxCharge"] = i_sector.second.bookCorrHistsX(
0728           "MaxCharge", "charge of strip with max. charge", "c_{cl,max}", "[APV counts]", 300, 75, -10., 290., "nph");
0729       i_sector.second.m_correlationHistsX["MaxIndex"] = i_sector.second.bookCorrHistsX(
0730           "MaxIndex", "cluster-index of strip with max. charge", "i_{cl,max}", "[# strips]", 10, 10, 0., 10., "nph");
0731       i_sector.second.m_correlationHistsX["ChargeOnEdges"] =
0732           i_sector.second.bookCorrHistsX("ChargeOnEdges",
0733                                          "fraction of charge on edge strips",
0734                                          "(c_{st,L}+c_{st,R})/c_{cl}",
0735                                          "",
0736                                          60,
0737                                          60,
0738                                          -0.1,
0739                                          1.1,
0740                                          "nph");
0741       i_sector.second.m_correlationHistsX["ChargeAsymmetry"] =
0742           i_sector.second.bookCorrHistsX("ChargeAsymmetry",
0743                                          "asymmetry of charge on edge strips",
0744                                          "(c_{st,L}-c_{st,R})/c_{cl}",
0745                                          "",
0746                                          110,
0747                                          55,
0748                                          -1.1,
0749                                          1.1,
0750                                          "nph");
0751       i_sector.second.m_correlationHistsX["ChargeLRplus"] =
0752           i_sector.second.bookCorrHistsX("ChargeLRplus",
0753                                          "fraction of charge not on maxStrip",
0754                                          "(c_{cl,L}+c_{cl,R})/c_{cl}",
0755                                          "",
0756                                          60,
0757                                          60,
0758                                          -0.1,
0759                                          1.1,
0760                                          "nph");
0761       i_sector.second.m_correlationHistsX["ChargeLRminus"] =
0762           i_sector.second.bookCorrHistsX("ChargeLRminus",
0763                                          "asymmetry of charge L and R of maxStrip",
0764                                          "(c_{cl,L}-c_{cl,R})/c_{cl}",
0765                                          "",
0766                                          110,
0767                                          55,
0768                                          -1.1,
0769                                          1.1,
0770                                          "nph");
0771       i_sector.second.m_correlationHistsX["SOverN"] =
0772           i_sector.second.bookCorrHistsX("SOverN", "signal over noise", "s/N", "", 100, 50, 0, sOverNMax, "nph");
0773       i_sector.second.m_correlationHistsX["WidthProj"] = i_sector.second.bookCorrHistsX(
0774           "WidthProj", "projected width", "w_{p}", "[# strips]", 200, 20, 0., widthMax, "nph");
0775       i_sector.second.m_correlationHistsX["WidthDiff"] = i_sector.second.bookCorrHistsX("WidthDiff",
0776                                                                                         "width difference",
0777                                                                                         "w_{p} - w_{cl}",
0778                                                                                         "[# strips]",
0779                                                                                         200,
0780                                                                                         20,
0781                                                                                         -widthMax / 2.,
0782                                                                                         widthMax / 2.,
0783                                                                                         "nph");
0784 
0785       i_sector.second.WidthVsWidthProjected = secDir.make<TH2F>("h2_widthVsWidthProj",
0786                                                                 "w_{cl} vs. w_{p};w_{p}  [# strips];w_{cl}  [# strips]",
0787                                                                 static_cast<int>(widthMax),
0788                                                                 0,
0789                                                                 widthMax,
0790                                                                 static_cast<int>(widthMax),
0791                                                                 0,
0792                                                                 widthMax);
0793       i_sector.second.PWidthVsWidthProjected =
0794           secDir.make<TProfile>("p_widthVsWidthProj",
0795                                 "w_{cl} vs. w_{p};w_{p}  [# strips];w_{cl}  [# strips]",
0796                                 static_cast<int>(widthMax),
0797                                 0,
0798                                 widthMax);
0799 
0800       i_sector.second.WidthDiffVsMaxStrip =
0801           secDir.make<TH2F>("h2_widthDiffVsMaxStrip",
0802                             "(w_{p} - w_{cl}) vs. n_{cl,max};n_{cl,max};w_{p} - w_{cl}  [# strips]",
0803                             800,
0804                             -10.,
0805                             790.,
0806                             static_cast<int>(widthMax),
0807                             -widthMax / 2.,
0808                             widthMax / 2.);
0809       i_sector.second.PWidthDiffVsMaxStrip =
0810           secDir.make<TProfile>("p_widthDiffVsMaxStrip",
0811                                 "(w_{p} - w_{cl}) vs. n_{cl,max};n_{cl,max};w_{p} - w_{cl}  [# strips]",
0812                                 800,
0813                                 -10.,
0814                                 790.);
0815 
0816       i_sector.second.WidthDiffVsSigmaXHit =
0817           secDir.make<TH2F>("h2_widthDiffVsSigmaXHit",
0818                             "(w_{p} - w_{cl}) vs. #sigma_{hit,x};#sigma_{hit,x}  [cm];w_{p} - w_{cl}  [# strips]",
0819                             100,
0820                             0.,
0821                             sigmaXMax,
0822                             100,
0823                             -10.,
0824                             10.);
0825       i_sector.second.PWidthDiffVsSigmaXHit =
0826           secDir.make<TProfile>("p_widthDiffVsSigmaXHit",
0827                                 "(w_{p} - w_{cl}) vs. #sigma_{hit,x};#sigma_{hit,x}  [cm];w_{p} - w_{cl}  [# strips]",
0828                                 100,
0829                                 0.,
0830                                 sigmaXMax);
0831 
0832       i_sector.second.WidthVsPhiSensX =
0833           secDir.make<TH2F>("h2_widthVsPhiSensX",
0834                             "w_{cl} vs. #phi_{module,x};#phi_{module,x}  [ ^{o}];w_{cl}  [# strips]",
0835                             93,
0836                             -93,
0837                             93,
0838                             static_cast<int>(widthMax),
0839                             0,
0840                             widthMax);
0841       i_sector.second.PWidthVsPhiSensX = secDir.make<TProfile>(
0842           "p_widthVsPhiSensX", "w_{cl} vs. #phi_{module,x};#phi_{module,x}  [ ^{o}];w_{cl}  [# strips]", 93, -93, 93);
0843     }
0844 
0845     // Hit Parameters (transform errors and residuals from [cm] in [mum])
0846     i_sector.second.m_correlationHistsX["SigmaXHit"] = i_sector.second.bookCorrHistsX(
0847         "SigmaXHit", "hit error", "#sigma_{hit,x}", "[#mum]", 105, 20, sigmaXMin * 10000., sigmaXMax * 10000., "np");
0848     i_sector.second.m_correlationHistsX["SigmaXTrk"] = i_sector.second.bookCorrHistsX(
0849         "SigmaXTrk", "track error", "#sigma_{trk,x}", "[#mum]", 105, 20, sigmaXMin * 10000., sigmaXMax * 10000., "np");
0850     i_sector.second.m_correlationHistsX["SigmaX"] = i_sector.second.bookCorrHistsX(
0851         "SigmaX", "residual error", "#sigma_{r,x}", "[#mum]", 105, 20, sigmaXMin * 10000., sigmaXMax * 10000., "np");
0852     i_sector.second.m_correlationHistsX["PhiSens"] = i_sector.second.bookCorrHistsX(
0853         "PhiSens", "track angle on sensor", "#phi_{module}", "[ ^{o}]", 96, 48, -3, 93, "nphtr");
0854     i_sector.second.m_correlationHistsX["PhiSensX"] = i_sector.second.bookCorrHistsX(
0855         "PhiSensX", "track angle on sensor", "#phi_{module,x}", "[ ^{o}]", 186, 93, -phiSensXMax, phiSensXMax, "nphtr");
0856     i_sector.second.m_correlationHistsX["PhiSensY"] = i_sector.second.bookCorrHistsX(
0857         "PhiSensY", "track angle on sensor", "#phi_{module,y}", "[ ^{o}]", 186, 93, -93, 93, "nphtr");
0858 
0859     i_sector.second.XHit = secDir.make<TH1F>("h_XHit", " hit measurement x_{hit};x_{hit}  [cm];# hits", 100, -20, 20);
0860     i_sector.second.XTrk = secDir.make<TH1F>("h_XTrk", "track prediction x_{trk};x_{trk}  [cm];# hits", 100, -20, 20);
0861     i_sector.second.SigmaX2 =
0862         secDir.make<TH1F>("h_SigmaX2",
0863                           "squared residual error #sigma_{r,x}^{2};#sigma_{r,x}^{2}  [#mum^{2}];# hits",
0864                           105,
0865                           sigmaXMin * 10000.,
0866                           sigmaX2Max * 10000. * 10000.);  //no mistake !
0867     i_sector.second.ResX = secDir.make<TH1F>(
0868         "h_ResX", "residual r_{x};x_{trk}-x_{hit}  [#mum];# hits", 100, -resXAbsMax * 10000., resXAbsMax * 10000.);
0869     i_sector.second.NorResX =
0870         secDir.make<TH1F>("h_NorResX",
0871                           "normalized residual r_{x}/#sigma_{r,x};(x_{trk}-x_{hit})/#sigma_{r,x};# hits",
0872                           100,
0873                           -norResXAbsMax,
0874                           norResXAbsMax);
0875     i_sector.second.ProbX = secDir.make<TH1F>(
0876         "h_ProbX", "residual probability;prob(r_{x}^{2}/#sigma_{r,x}^{2},1);# hits", 60, probXMin, probXMax);
0877 
0878     i_sector.second.PhiSensXVsBarycentreX =
0879         secDir.make<TH2F>("h2_phiSensXVsBarycentreX",
0880                           "#phi_{module,x} vs. b_{cl,x};b_{cl,x}  [# channels];#phi_{module,x}  [ ^{o}]",
0881                           200,
0882                           -10.,
0883                           790.,
0884                           93,
0885                           -93,
0886                           93);
0887     i_sector.second.PPhiSensXVsBarycentreX =
0888         secDir.make<TProfile>("p_phiSensXVsBarycentreX",
0889                               "#phi_{module,x} vs. b_{cl,x};b_{cl,x}  [# channels];#phi_{module,x}  [ ^{o}]",
0890                               200,
0891                               -10.,
0892                               790.);
0893 
0894     if (pixelSector) {
0895       i_sector.second.m_correlationHistsY["SigmaYHit"] = i_sector.second.bookCorrHistsY(
0896           "SigmaYHit", "hit error", "#sigma_{hit,y}", "[#mum]", 105, 20, sigmaXMin * 10000., sigmaXMax * 10000., "np");
0897       i_sector.second.m_correlationHistsY["SigmaYTrk"] = i_sector.second.bookCorrHistsY(
0898           "SigmaYTrk", "track error", "#sigma_{trk,y}", "[#mum]", 105, 20, sigmaXMin * 10000., sigmaXMax * 10000., "np");
0899       i_sector.second.m_correlationHistsY["SigmaY"] = i_sector.second.bookCorrHistsY(
0900           "SigmaY", "residual error", "#sigma_{r,y}", "[#mum]", 105, 20, sigmaXMin * 10000., sigmaXMax * 10000., "np");
0901       i_sector.second.m_correlationHistsY["PhiSens"] = i_sector.second.bookCorrHistsY(
0902           "PhiSens", "track angle on sensor", "#phi_{module}", "[ ^{o}]", 96, 48, -3, 93, "nphtr");
0903       i_sector.second.m_correlationHistsY["PhiSensX"] = i_sector.second.bookCorrHistsY("PhiSensX",
0904                                                                                        "track angle on sensor",
0905                                                                                        "#phi_{module,x}",
0906                                                                                        "[ ^{o}]",
0907                                                                                        186,
0908                                                                                        93,
0909                                                                                        -phiSensXMax,
0910                                                                                        phiSensXMax,
0911                                                                                        "nphtr");
0912       i_sector.second.m_correlationHistsY["PhiSensY"] = i_sector.second.bookCorrHistsY(
0913           "PhiSensY", "track angle on sensor", "#phi_{module,y}", "[ ^{o}]", 186, 93, -93, 93, "nphtr");
0914 
0915       i_sector.second.YHit = secDir.make<TH1F>("h_YHit", " hit measurement y_{hit};y_{hit}  [cm];# hits", 100, -20, 20);
0916       i_sector.second.YTrk = secDir.make<TH1F>("h_YTrk", "track prediction y_{trk};y_{trk}  [cm];# hits", 100, -20, 20);
0917       i_sector.second.SigmaY2 =
0918           secDir.make<TH1F>("h_SigmaY2",
0919                             "squared residual error #sigma_{r,y}^{2};#sigma_{r,y}^{2}  [#mum^{2}];# hits",
0920                             105,
0921                             sigmaXMin * 10000.,
0922                             sigmaX2Max * 10000. * 10000.);  //no mistake !
0923       i_sector.second.ResY = secDir.make<TH1F>(
0924           "h_ResY", "residual r_{y};y_{trk}-y_{hit}  [#mum];# hits", 100, -resXAbsMax * 10000., resXAbsMax * 10000.);
0925       i_sector.second.NorResY =
0926           secDir.make<TH1F>("h_NorResY",
0927                             "normalized residual r_{y}/#sigma_{r,y};(y_{trk}-y_{hit})/#sigma_{r,y};# hits",
0928                             100,
0929                             -norResXAbsMax,
0930                             norResXAbsMax);
0931       i_sector.second.ProbY = secDir.make<TH1F>(
0932           "h_ProbY", "residual probability;prob(r_{y}^{2}/#sigma_{r,y}^{2},1);# hits", 60, probXMin, probXMax);
0933 
0934       i_sector.second.PhiSensYVsBarycentreY =
0935           secDir.make<TH2F>("h2_phiSensYVsBarycentreY",
0936                             "#phi_{module,y} vs. b_{cl,y};b_{cl,y}  [# channels];#phi_{module,y}  [ ^{o}]",
0937                             200,
0938                             -10.,
0939                             790.,
0940                             93,
0941                             -93,
0942                             93);
0943       i_sector.second.PPhiSensYVsBarycentreY =
0944           secDir.make<TProfile>("p_phiSensYVsBarycentreY",
0945                                 "#phi_{module,y} vs. b_{cl,y};b_{cl,y}  [# channels];#phi_{module,y}  [ ^{o}]",
0946                                 200,
0947                                 -10.,
0948                                 790.);
0949     }
0950 
0951     // Track Parameters
0952     i_sector.second.m_correlationHistsX["HitsValid"] =
0953         i_sector.second.bookCorrHistsX("HitsValid", "# hits", "[valid]", 50, 0, 50, "npt");
0954     i_sector.second.m_correlationHistsX["HitsInvalid"] =
0955         i_sector.second.bookCorrHistsX("HitsInvalid", "# hits", "[invalid]", 20, 0, 20, "npt");
0956     i_sector.second.m_correlationHistsX["Hits2D"] =
0957         i_sector.second.bookCorrHistsX("Hits2D", "# hits", "[2D]", 20, 0, 20, "npt");
0958     i_sector.second.m_correlationHistsX["LayersMissed"] =
0959         i_sector.second.bookCorrHistsX("LayersMissed", "# layers", "[missed]", 10, 0, 10, "npt");
0960     i_sector.second.m_correlationHistsX["HitsPixel"] =
0961         i_sector.second.bookCorrHistsX("HitsPixel", "# hits", "[pixel]", 10, 0, 10, "npt");
0962     i_sector.second.m_correlationHistsX["HitsStrip"] =
0963         i_sector.second.bookCorrHistsX("HitsStrip", "# hits", "[strip]", 40, 0, 40, "npt");
0964     i_sector.second.m_correlationHistsX["HitsGood"] =
0965         i_sector.second.bookCorrHistsX("HitsGood", "# hits", "[good]", 50, 0, 50, "npt");
0966     i_sector.second.m_correlationHistsX["NorChi2"] =
0967         i_sector.second.bookCorrHistsX("NorChi2", "#chi^{2}/f", "", 50, 0, norChi2Max, "npr");
0968     i_sector.second.m_correlationHistsX["Theta"] =
0969         i_sector.second.bookCorrHistsX("Theta", "#theta", "[ ^{o}]", 40, -10, 190, "npt");
0970     i_sector.second.m_correlationHistsX["Phi"] =
0971         i_sector.second.bookCorrHistsX("Phi", "#phi", "[ ^{o}]", 76, -190, 190, "npt");
0972     i_sector.second.m_correlationHistsX["D0Beamspot"] =
0973         i_sector.second.bookCorrHistsX("D0Beamspot", "d_{0, BS}", "[cm]", 40, -d0Max, d0Max, "npt");
0974     i_sector.second.m_correlationHistsX["Dz"] =
0975         i_sector.second.bookCorrHistsX("Dz", "d_{z}", "[cm]", 40, -dzMax, dzMax, "npt");
0976     i_sector.second.m_correlationHistsX["Pt"] =
0977         i_sector.second.bookCorrHistsX("Pt", "p_{t}", "[GeV]", 50, 0, pMax, "npt");
0978     i_sector.second.m_correlationHistsX["P"] = i_sector.second.bookCorrHistsX("P", "|p|", "[GeV]", 50, 0, pMax, "npt");
0979     i_sector.second.m_correlationHistsX["InvP"] =
0980         i_sector.second.bookCorrHistsX("InvP", "1/|p|", "[GeV^{-1}]", 25, 0, invPMax, "t");
0981     i_sector.second.m_correlationHistsX["MeanAngle"] =
0982         i_sector.second.bookCorrHistsX("MeanAngle", "<#phi_{module}>", "[ ^{o}]", 25, -5, 95, "npt");
0983     //i_sector.second.m_correlationHistsX[""] = i_sector.second.bookCorrHistsX("","","",,,,"nphtr");
0984 
0985     if (pixelSector) {
0986       i_sector.second.m_correlationHistsY["HitsValid"] =
0987           i_sector.second.bookCorrHistsY("HitsValid", "# hits", "[valid]", 50, 0, 50, "npt");
0988       i_sector.second.m_correlationHistsY["HitsInvalid"] =
0989           i_sector.second.bookCorrHistsY("HitsInvalid", "# hits", "[invalid]", 20, 0, 20, "npt");
0990       i_sector.second.m_correlationHistsY["Hits2D"] =
0991           i_sector.second.bookCorrHistsY("Hits2D", "# hits", "[2D]", 20, 0, 20, "npt");
0992       i_sector.second.m_correlationHistsY["LayersMissed"] =
0993           i_sector.second.bookCorrHistsY("LayersMissed", "# layers", "[missed]", 10, 0, 10, "npt");
0994       i_sector.second.m_correlationHistsY["HitsPixel"] =
0995           i_sector.second.bookCorrHistsY("HitsPixel", "# hits", "[pixel]", 10, 0, 10, "npt");
0996       i_sector.second.m_correlationHistsY["HitsStrip"] =
0997           i_sector.second.bookCorrHistsY("HitsStrip", "# hits", "[strip]", 40, 0, 40, "npt");
0998       i_sector.second.m_correlationHistsY["HitsGood"] =
0999           i_sector.second.bookCorrHistsY("HitsGood", "# hits", "[good]", 50, 0, 50, "npt");
1000       i_sector.second.m_correlationHistsY["NorChi2"] =
1001           i_sector.second.bookCorrHistsY("NorChi2", "#chi^{2}/f", "", 50, 0, norChi2Max, "npr");
1002       i_sector.second.m_correlationHistsY["Theta"] =
1003           i_sector.second.bookCorrHistsY("Theta", "#theta", "[ ^{o}]", 40, -10, 190, "npt");
1004       i_sector.second.m_correlationHistsY["Phi"] =
1005           i_sector.second.bookCorrHistsY("Phi", "#phi", "[ ^{o}]", 76, -190, 190, "npt");
1006       i_sector.second.m_correlationHistsY["D0Beamspot"] =
1007           i_sector.second.bookCorrHistsY("D0Beamspot", "d_{0, BS}", "[cm]", 40, -d0Max, d0Max, "npt");
1008       i_sector.second.m_correlationHistsY["Dz"] =
1009           i_sector.second.bookCorrHistsY("Dz", "d_{z}", "[cm]", 40, -dzMax, dzMax, "npt");
1010       i_sector.second.m_correlationHistsY["Pt"] =
1011           i_sector.second.bookCorrHistsY("Pt", "p_{t}", "[GeV]", 50, 0, pMax, "npt");
1012       i_sector.second.m_correlationHistsY["P"] =
1013           i_sector.second.bookCorrHistsY("P", "|p|", "[GeV]", 50, 0, pMax, "npt");
1014       i_sector.second.m_correlationHistsY["InvP"] =
1015           i_sector.second.bookCorrHistsY("InvP", "1/|p|", "[GeV^{-1}]", 25, 0, invPMax, "t");
1016       i_sector.second.m_correlationHistsY["MeanAngle"] =
1017           i_sector.second.bookCorrHistsY("MeanAngle", "<#phi_{module}>", "[ ^{o}]", 25, -5, 95, "npt");
1018     }
1019 
1020     // (transform errors and residuals from [cm] in [mum])
1021     for (auto const& i_errHists : v_errHists) {
1022       double xMin(0.01 * (i_errHists - 1)), xMax(0.01 * (i_errHists));
1023       std::stringstream sigmaXHit, sigmaXTrk, sigmaX;
1024       sigmaXHit << "h_sigmaXHit_" << i_errHists;
1025       sigmaXTrk << "h_sigmaXTrk_" << i_errHists;
1026       sigmaX << "h_sigmaX_" << i_errHists;
1027       i_sector.second.m_sigmaX["sigmaXHit"].push_back(
1028           secDir.make<TH1F>(sigmaXHit.str().c_str(),
1029                             "hit error #sigma_{hit,x};#sigma_{hit,x}  [#mum];# hits",
1030                             100,
1031                             xMin * 10000.,
1032                             xMax * 10000.));
1033       i_sector.second.m_sigmaX["sigmaXTrk"].push_back(
1034           secDir.make<TH1F>(sigmaXTrk.str().c_str(),
1035                             "track error #sigma_{trk,x};#sigma_{trk,x}  [#mum];# hits",
1036                             100,
1037                             xMin * 10000.,
1038                             xMax * 10000.));
1039       i_sector.second.m_sigmaX["sigmaX"].push_back(
1040           secDir.make<TH1F>(sigmaX.str().c_str(),
1041                             "residual error #sigma_{r,x};#sigma_{r,x}  [#mum];# hits",
1042                             100,
1043                             xMin * 10000.,
1044                             xMax * 10000.));
1045       if (pixelSector) {
1046         std::stringstream sigmaYHit, sigmaYTrk, sigmaY;
1047         sigmaYHit << "h_sigmaYHit_" << i_errHists;
1048         sigmaYTrk << "h_sigmaYTrk_" << i_errHists;
1049         sigmaY << "h_sigmaY_" << i_errHists;
1050         i_sector.second.m_sigmaY["sigmaYHit"].push_back(
1051             secDir.make<TH1F>(sigmaYHit.str().c_str(),
1052                               "hit error #sigma_{hit,y};#sigma_{hit,y}  [#mum];# hits",
1053                               100,
1054                               xMin * 10000.,
1055                               xMax * 10000.));
1056         i_sector.second.m_sigmaY["sigmaYTrk"].push_back(
1057             secDir.make<TH1F>(sigmaYTrk.str().c_str(),
1058                               "track error #sigma_{trk,y};#sigma_{trk,y}  [#mum];# hits",
1059                               100,
1060                               xMin * 10000.,
1061                               xMax * 10000.));
1062         i_sector.second.m_sigmaY["sigmaY"].push_back(
1063             secDir.make<TH1F>(sigmaY.str().c_str(),
1064                               "residual error #sigma_{r,y};#sigma_{r,y}  [#mum];# hits",
1065                               100,
1066                               xMin * 10000.,
1067                               xMax * 10000.));
1068       }
1069     }
1070   }
1071 }
1072 
1073 void ApeEstimator::bookSectorHistsForApeCalculation() {
1074   std::vector<unsigned int> v_errHists(parameterSet_.getParameter<std::vector<unsigned int> >("vErrHists"));
1075   for (std::vector<unsigned int>::iterator i_errHists = v_errHists.begin(); i_errHists != v_errHists.end();
1076        ++i_errHists) {
1077     for (std::vector<unsigned int>::iterator i_errHists2 = i_errHists; i_errHists2 != v_errHists.end();) {
1078       ++i_errHists2;
1079       if (*i_errHists == *i_errHists2) {
1080         edm::LogError("BookSectorHists") << "Value of vErrHists in config exists twice: " << *i_errHists
1081                                          << "\n... delete one of both";
1082         v_errHists.erase(i_errHists2);
1083       }
1084     }
1085   }
1086 
1087   for (auto& i_sector : m_tkSector_) {
1088     edm::Service<TFileService> fileService;
1089     if (!fileService) {
1090       throw edm::Exception(edm::errors::Configuration, "TFileService is not registered in cfg file");
1091     }
1092 
1093     std::stringstream sector;
1094     sector << "Sector_" << i_sector.first;
1095     TFileDirectory secDir = fileService->mkdir(sector.str());
1096 
1097     // Dummy histo containing the sector name as title
1098     i_sector.second.Name = secDir.make<TH1F>("z_name", i_sector.second.name.c_str(), 1, 0, 1);
1099 
1100     // Do not book histos for empty sectors
1101     if (i_sector.second.v_rawId.empty()) {
1102       continue;
1103     }
1104 
1105     // Distributions in each interval (stay in [cm], to have all calculations in [cm])
1106     if (m_resErrBins_
1107             .empty()) {  // default if no selection taken into account: calculate APE with one bin with residual error 0-100um
1108       m_resErrBins_[1].first = 0.;
1109       m_resErrBins_[1].second = 0.01;
1110     }
1111     for (auto const& i_errBins : m_resErrBins_) {
1112       std::stringstream interval;
1113       interval << "Interval_" << i_errBins.first;
1114       TFileDirectory intDir = secDir.mkdir(interval.str());
1115       i_sector.second.m_binnedHists[i_errBins.first]["sigmaX"] =
1116           intDir.make<TH1F>("h_sigmaX", "residual resolution #sigma_{x};#sigma_{x}  [cm];# hits", 100, 0., 0.01);
1117       i_sector.second.m_binnedHists[i_errBins.first]["norResX"] = intDir.make<TH1F>(
1118           "h_norResX", "normalized residual r_{x}/#sigma_{r,x};(x_{trk}-x_{hit})/#sigma_{r,x};# hits", 100, -10, 10);
1119       if (i_sector.second.isPixel) {
1120         i_sector.second.m_binnedHists[i_errBins.first]["sigmaY"] =
1121             intDir.make<TH1F>("h_sigmaY", "residual resolution #sigma_{y};#sigma_{y}  [cm];# hits", 100, 0., 0.01);
1122         i_sector.second.m_binnedHists[i_errBins.first]["norResY"] = intDir.make<TH1F>(
1123             "h_norResY", "normalized residual r_{y}/#sigma_{r,y};(y_{trk}-y_{hit})/#sigma_{r,y};# hits", 100, -10, 10);
1124       }
1125     }
1126 
1127     TFileDirectory resDir = secDir.mkdir("Results");
1128 
1129     // TTree containing rawIds of all modules in sector
1130     unsigned int rawId(0);
1131     i_sector.second.RawId = resDir.make<TTree>("rawIdTree", "Tree containing rawIds of all modules in sector");
1132     i_sector.second.RawId->Branch("RawId", &rawId, "RawId/i");
1133     for (auto const& i_rawId : i_sector.second.v_rawId) {
1134       rawId = i_rawId;
1135       i_sector.second.RawId->Fill();
1136     }
1137 
1138     // Result plots (one hist per sector containing one bin per interval)
1139     // (transform errors and residuals from [cm] in [mum])
1140     std::vector<double> v_binX(parameterSet_.getParameter<std::vector<double> >("residualErrorBinning"));
1141     for (auto& i_binX : v_binX) {
1142       i_binX *= 10000.;
1143     }
1144     i_sector.second.EntriesX =
1145         resDir.make<TH1F>("h_entriesX", "# hits used;#sigma_{x}  [#mum];# hits", v_binX.size() - 1, &(v_binX[0]));
1146     if (i_sector.second.isPixel) {
1147       i_sector.second.EntriesY =
1148           resDir.make<TH1F>("h_entriesY", "# hits used;#sigma_{y}  [#mum];# hits", v_binX.size() - 1, &(v_binX[0]));
1149     }
1150 
1151     // In fact these are un-needed Analyzer plots, but I want to have them always for every sector visible
1152     // (transform errors and residuals from [cm] in [mum])
1153     i_sector.second.ResX = resDir.make<TH1F>(
1154         "h_ResX", "residual r_{x};x_{trk}-x_{hit}  [#mum];# hits", 100, -0.03 * 10000., 0.03 * 10000.);
1155     i_sector.second.NorResX = resDir.make<TH1F>(
1156         "h_NorResX", "normalized residual r_{x}/#sigma_{r,x};(x_{trk}-x_{hit})/#sigma_{r,x};# hits", 100, -5., 5.);
1157     if (i_sector.second.isPixel) {
1158       i_sector.second.ResY = resDir.make<TH1F>(
1159           "h_ResY", "residual r_{y};y_{trk}-y_{hit}  [#mum];# hits", 100, -0.03 * 10000., 0.03 * 10000.);
1160       i_sector.second.NorResY = resDir.make<TH1F>(
1161           "h_NorResY", "normalized residual r_{y}/#sigma_{r,y};(y_{trk}-y_{hit})/#sigma_{r,y};# hits", 100, -5., 5.);
1162     }
1163   }
1164 }
1165 
1166 // -----------------------------------------------------------------------------------------------------------
1167 
1168 void ApeEstimator::bookTrackHists() {
1169   bool zoomHists(parameterSet_.getParameter<bool>("zoomHists"));
1170 
1171   int trackSizeBins = zoomHists ? 6 : 201;
1172   double trackSizeMax = trackSizeBins - 1;
1173 
1174   double chi2Max = zoomHists ? 100. : 2000.;
1175   double norChi2Max = zoomHists ? 5. : 1000.;
1176   double d0max = zoomHists ? 0.02 : 40.;  // cosmics: 100.|100.
1177   double dzmax = zoomHists ? 15. : 100.;  // cosmics: 200.|600.
1178   double pMax = zoomHists ? 200. : 2000.;
1179 
1180   edm::Service<TFileService> fileService;
1181   TFileDirectory evtDir = fileService->mkdir("EventVariables");
1182   tkDetector_.TrkSize =
1183       evtDir.make<TH1F>("h_trackSize", "# tracks  [all];# tracks;# events", trackSizeBins, -1, trackSizeMax);
1184   tkDetector_.TrkSizeGood =
1185       evtDir.make<TH1F>("h_trackSizeGood", "# tracks  [good];# tracks;# events", trackSizeBins, -1, trackSizeMax);
1186   TFileDirectory trkDir = fileService->mkdir("TrackVariables");
1187   tkDetector_.HitsSize = trkDir.make<TH1F>("h_hitsSize", "# hits;# hits;# tracks", 51, -1, 50);
1188   tkDetector_.HitsValid = trkDir.make<TH1F>("h_hitsValid", "# hits  [valid];# hits  [valid];# tracks", 51, -1, 50);
1189   tkDetector_.HitsInvalid =
1190       trkDir.make<TH1F>("h_hitsInvalid", "# hits  [invalid];# hits  [invalid];# tracks", 21, -1, 20);
1191   tkDetector_.Hits2D = trkDir.make<TH1F>("h_hits2D", "# hits  [2D];# hits  [2D];# tracks", 21, -1, 20);
1192   tkDetector_.LayersMissed =
1193       trkDir.make<TH1F>("h_layersMissed", "# layers  [missed];# layers  [missed];# tracks", 11, -1, 10);
1194   tkDetector_.HitsPixel = trkDir.make<TH1F>("h_hitsPixel", "# hits  [pixel];# hits  [pixel];# tracks", 11, -1, 10);
1195   tkDetector_.HitsStrip = trkDir.make<TH1F>("h_hitsStrip", "# hits  [strip];# hits  [strip];# tracks", 41, -1, 40);
1196   tkDetector_.Charge = trkDir.make<TH1F>("h_charge", "charge q;q  [e];# tracks", 5, -2, 3);
1197   tkDetector_.Chi2 = trkDir.make<TH1F>("h_chi2", " #chi^{2};#chi^{2};# tracks", 100, 0, chi2Max);
1198   tkDetector_.Ndof = trkDir.make<TH1F>("h_ndof", "# degrees of freedom f;f;# tracks", 101, -1, 100);
1199   tkDetector_.NorChi2 = trkDir.make<TH1F>("h_norChi2", "normalized #chi^{2};#chi^{2}/f;# tracks", 200, 0, norChi2Max);
1200   tkDetector_.Prob = trkDir.make<TH1F>("h_prob", " #chi^{2} probability;prob(#chi^{2},f);# tracks", 50, 0, 1);
1201   tkDetector_.Eta = trkDir.make<TH1F>("h_eta", "pseudorapidity #eta;#eta;# tracks", 100, -5, 5);
1202   tkDetector_.EtaErr = trkDir.make<TH1F>("h_etaErr", "Error of #eta;#sigma(#eta);# tracks", 100, 0, 0.001);
1203   tkDetector_.EtaSig =
1204       trkDir.make<TH1F>("h_etaSig", "Significance of #eta;#eta/#sigma(#eta);# tracks", 100, -20000, 20000);
1205   tkDetector_.Theta = trkDir.make<TH1F>("h_theta", "polar angle #theta;#theta  [ ^{o}];# tracks", 100, -10, 190);
1206   tkDetector_.Phi = trkDir.make<TH1F>("h_phi", "azimuth angle #phi;#phi  [ ^{o}];# tracks", 190, -190, 190);
1207   tkDetector_.PhiErr = trkDir.make<TH1F>("h_phiErr", "Error of #phi;#sigma(#phi)  [ ^{o}];# tracks", 100, 0, 0.04);
1208   tkDetector_.PhiSig =
1209       trkDir.make<TH1F>("h_phiSig", "Significance of #phi;#phi/#sigma(#phi)  [ ^{o}];# tracks", 100, -50000, 50000);
1210   tkDetector_.D0Beamspot = trkDir.make<TH1F>(
1211       "h_d0Beamspot", "Closest approach d_{0} wrt. beamspot;d_{0, BS}  [cm];# tracks", 200, -d0max, d0max);
1212   tkDetector_.D0BeamspotErr =
1213       trkDir.make<TH1F>("h_d0BeamspotErr", "Error of d_{0, BS};#sigma(d_{0, BS})  [cm];# tracks", 200, 0, 0.01);
1214   tkDetector_.D0BeamspotSig = trkDir.make<TH1F>(
1215       "h_d0BeamspotSig", "Significance of d_{0, BS};d_{0, BS}/#sigma(d_{0, BS});# tracks", 100, -5, 5);
1216   tkDetector_.Dz = trkDir.make<TH1F>("h_dz", "Closest approach d_{z};d_{z}  [cm];# tracks", 200, -dzmax, dzmax);
1217   tkDetector_.DzErr = trkDir.make<TH1F>("h_dzErr", "Error of d_{z};#sigma(d_{z})  [cm];# tracks", 200, 0, 0.01);
1218   tkDetector_.DzSig =
1219       trkDir.make<TH1F>("h_dzSig", "Significance of d_{z};d_{z}/#sigma(d_{z});# tracks", 100, -10000, 10000);
1220   tkDetector_.Pt = trkDir.make<TH1F>("h_pt", "transverse momentum p_{t};p_{t}  [GeV];# tracks", 100, 0, pMax);
1221   tkDetector_.PtErr = trkDir.make<TH1F>("h_ptErr", "Error of p_{t};#sigma(p_{t})  [GeV];# tracks", 100, 0, 1.6);
1222   tkDetector_.PtSig = trkDir.make<TH1F>("h_ptSig", "Significance of p_{t};p_{t}/#sigma(p_{t});# tracks", 100, 0, 200);
1223   tkDetector_.P = trkDir.make<TH1F>("h_p", "momentum magnitude |p|;|p|  [GeV];# tracks", 100, 0, pMax);
1224   tkDetector_.MeanAngle = trkDir.make<TH1F>(
1225       "h_meanAngle", "mean angle on module <#phi_{module}>;<#phi_{module}>  [ ^{o}];# tracks", 100, -5, 95);
1226   tkDetector_.HitsGood = trkDir.make<TH1F>("h_hitsGood", "# hits  [good];# hits  [good];# tracks", 51, -1, 50);
1227 
1228   tkDetector_.MeanAngleVsHits = trkDir.make<TH2F>(
1229       "h2_meanAngleVsHits", "<#phi_{module}> vs. # hits;# hits;<#phi_{module}>  [ ^{o}]", 51, -1, 50, 50, -5, 95);
1230   tkDetector_.HitsGoodVsHitsValid =
1231       trkDir.make<TH2F>("h2_hitsGoodVsHitsValid",
1232                         "# hits  [good] vs. # hits  [valid];# hits  [valid];# hits  [good]",
1233                         51,
1234                         -1,
1235                         50,
1236                         51,
1237                         -1,
1238                         50);
1239   tkDetector_.HitsPixelVsEta =
1240       trkDir.make<TH2F>("h2_hitsPixelVsEta", "# hits  [pixel] vs. #eta;#eta;# hits  [pixel]", 60, -3, 3, 11, -1, 10);
1241   tkDetector_.HitsPixelVsTheta = trkDir.make<TH2F>(
1242       "h2_hitsPixelVsTheta", "# hits  [pixel] vs. #theta;#theta;# hits  [pixel]", 100, -10, 190, 11, -1, 10);
1243   tkDetector_.HitsStripVsEta =
1244       trkDir.make<TH2F>("h2_hitsStripVsEta", "# hits  [strip] vs. #eta;#eta;# hits  [strip]", 60, -3, 3, 31, -1, 40);
1245   tkDetector_.HitsStripVsTheta = trkDir.make<TH2F>(
1246       "h2_hitsStripVsTheta", "# hits  [strip] vs. #theta;#theta;# hits  [strip]", 100, -10, 190, 31, -1, 40);
1247   tkDetector_.PtVsEta = trkDir.make<TH2F>("h2_ptVsEta", "p_{t} vs. #eta;#eta;p_{t}  [GeV]", 60, -3, 3, 100, 0, pMax);
1248   tkDetector_.PtVsTheta =
1249       trkDir.make<TH2F>("h2_ptVsTheta", "p_{t} vs. #theta;#theta;p_{t}  [GeV]", 100, -10, 190, 100, 0, pMax);
1250 
1251   tkDetector_.PMeanAngleVsHits = trkDir.make<TProfile>(
1252       "p_meanAngleVsHits", "<#phi_{module}> vs. # hits;# hits;<#phi_{module}>  [ ^{o}]", 51, -1, 50);
1253   tkDetector_.PHitsGoodVsHitsValid = trkDir.make<TProfile>(
1254       "p_hitsGoodVsHitsValid", "# hits  [good] vs. # hits  [valid];# hits  [valid];# hits  [good]", 51, -1, 50);
1255   tkDetector_.PHitsPixelVsEta =
1256       trkDir.make<TProfile>("p_hitsPixelVsEta", "# hits  [pixel] vs. #eta;#eta;# hits  [pixel]", 60, -3, 3);
1257   tkDetector_.PHitsPixelVsTheta =
1258       trkDir.make<TProfile>("p_hitsPixelVsTheta", "# hits  [pixel] vs. #theta;#theta;# hits  [pixel]", 100, -10, 190);
1259   tkDetector_.PHitsStripVsEta =
1260       trkDir.make<TProfile>("p_hitsStripVsEta", "# hits  [strip] vs. #eta;#eta;# hits  [strip]", 60, -3, 3);
1261   tkDetector_.PHitsStripVsTheta =
1262       trkDir.make<TProfile>("p_hitsStripVsTheta", "# hits  [strip] vs. #theta;#theta;# hits  [strip]", 100, -10, 190);
1263   tkDetector_.PPtVsEta = trkDir.make<TProfile>("p_ptVsEta", "p_{t} vs. #eta;#eta;p_{t}  [GeV]", 60, -3, 3);
1264   tkDetector_.PPtVsTheta = trkDir.make<TProfile>("p_ptVsTheta", "p_{t} vs. #theta;#theta;p_{t}  [GeV]", 100, -10, 190);
1265 }
1266 
1267 // -----------------------------------------------------------------------------------------------------------
1268 
1269 TrackStruct::TrackParameterStruct ApeEstimator::fillTrackVariables(const reco::Track& track,
1270                                                                    const Trajectory& traj,
1271                                                                    const reco::BeamSpot& beamSpot) {
1272   const math::XYZPoint beamPoint(beamSpot.x0(), beamSpot.y0(), beamSpot.z0());
1273   double d0BeamspotErr =
1274       std::sqrt(track.d0Error() * track.d0Error() + 0.5 * beamSpot.BeamWidthX() * beamSpot.BeamWidthX() +
1275                 0.5 * beamSpot.BeamWidthY() * beamSpot.BeamWidthY());
1276 
1277   const static TrajectoryStateCombiner tsoscomb;
1278 
1279   const reco::HitPattern& hitPattern(track.hitPattern());
1280 
1281   TrackStruct::TrackParameterStruct trkParams;
1282 
1283   trkParams.hitsSize = track.recHitsSize();
1284   trkParams.hitsValid = track.found();  // invalid is every hit from every single module that expects a hit
1285   trkParams.hitsInvalid = trkParams.hitsSize - trkParams.hitsValid;
1286   trkParams.layersMissed =
1287       track.lost();  // lost hit means, that a crossed layer doesn't contain a hit (can be more than one invalid hit)
1288   trkParams.hitsPixel = hitPattern.numberOfValidPixelHits();
1289   trkParams.hitsStrip = hitPattern.numberOfValidStripHits();
1290   trkParams.charge = track.charge();
1291   trkParams.chi2 = track.chi2();
1292   trkParams.ndof = track.ndof();
1293   trkParams.norChi2 = trkParams.chi2 / trkParams.ndof;
1294   trkParams.prob = TMath::Prob(trkParams.chi2, trkParams.ndof);
1295   trkParams.eta = track.eta();
1296   trkParams.etaErr = track.etaError();
1297   trkParams.theta = track.theta();
1298   trkParams.phi = track.phi();
1299   trkParams.phiErr = track.phiError();
1300   trkParams.d0 = track.d0();
1301   trkParams.d0Beamspot = -1. * track.dxy(beamPoint);
1302   trkParams.d0BeamspotErr = d0BeamspotErr;
1303   trkParams.dz = track.dz();
1304   trkParams.dzErr = track.dzError();
1305   trkParams.dzBeamspot = track.dz(beamPoint);
1306   trkParams.p = track.p();
1307   trkParams.pt = track.pt();
1308   trkParams.ptErr = track.ptError();
1309 
1310   const std::vector<TrajectoryMeasurement>& v_meas = traj.measurements();
1311 
1312   int count2D(0);
1313   float meanPhiSensToNorm(0.F);
1314   for (auto const& i_meas : v_meas) {
1315     const TrajectoryMeasurement& meas = i_meas;
1316     const TransientTrackingRecHit& hit = *meas.recHit();
1317     const TrackingRecHit& recHit = *hit.hit();
1318     if (this->isHit2D(recHit))
1319       ++count2D;
1320 
1321     TrajectoryStateOnSurface tsos = tsoscomb(meas.forwardPredictedState(), meas.backwardPredictedState());
1322     const align::LocalVector mom(tsos.localDirection());
1323     meanPhiSensToNorm += atan(fabs(sqrt(mom.x() * mom.x() + mom.y() * mom.y()) / mom.z()));
1324   }
1325   meanPhiSensToNorm *= (1. / static_cast<float>(trkParams.hitsSize));
1326 
1327   trkParams.hits2D = count2D;
1328   trkParams.meanPhiSensToNorm = meanPhiSensToNorm;
1329 
1330   if (parameterSet_.getParameter<bool>("applyTrackCuts")) {
1331     trackCut_ = false;
1332     if (trkParams.hitsStrip < 11 || trkParams.hits2D < 2 || trkParams.hitsPixel < 2 ||  //trkParams.hitsInvalid>2 ||
1333         trkParams.hitsStrip > 35 || trkParams.hitsPixel > 7 || trkParams.norChi2 > 5. || trkParams.pt < 25. ||
1334         trkParams.pt > 150. || std::abs(trkParams.d0Beamspot) > 0.02 || std::abs(trkParams.dz) > 15.)
1335       trackCut_ = true;
1336   } else {
1337     trackCut_ = false;
1338   }
1339 
1340   return trkParams;
1341 }
1342 
1343 TrackStruct::HitParameterStruct ApeEstimator::fillHitVariables(const TrajectoryMeasurement& i_meas,
1344                                                                const edm::EventSetup& iSetup) {
1345   TrackStruct::HitParameterStruct hitParams;
1346 
1347   const static TrajectoryStateCombiner tsoscomb;
1348 
1349   const TrajectoryMeasurement& meas = i_meas;
1350   const TransientTrackingRecHit& hit = *meas.recHit();
1351   const TrackingRecHit& recHit = *hit.hit();
1352   const TrajectoryStateOnSurface& tsos = tsoscomb(meas.forwardPredictedState(), meas.backwardPredictedState());
1353 
1354   const DetId& detId(hit.geographicalId());
1355   const DetId::Detector& detector = detId.det();
1356   if (detector != DetId::Tracker) {
1357     hitParams.hitState = TrackStruct::notInTracker;
1358     return hitParams;
1359   }
1360   const uint32_t rawId(detId.rawId());
1361 
1362   for (auto& i_sector : m_tkSector_) {
1363     for (auto const& i_rawId : i_sector.second.v_rawId) {
1364       if (rawId == i_rawId) {
1365         hitParams.v_sector.push_back(i_sector.first);
1366         break;
1367       }
1368     }
1369   }
1370 
1371   const align::LocalVector& mom(tsos.localDirection());
1372   int xMomentum(0), yMomentum(0), zMomentum(0);
1373   xMomentum = mom.x() > 0. ? 1 : -1;
1374   yMomentum = mom.y() > 0. ? 1 : -1;
1375   zMomentum = mom.z() > 0. ? 1 : -1;
1376   float phiSensX =
1377       std::atan(std::fabs(mom.x() / mom.z())) *
1378       static_cast<float>(
1379           m_tkTreeVar_[rawId].vDirection);  // check for orientation of E- and B- Field (thoughts for barrel)
1380   float phiSensY = std::atan(std::fabs(mom.y() / mom.z())) * static_cast<float>(m_tkTreeVar_[rawId].vDirection);
1381   hitParams.phiSens = std::atan(std::fabs(std::sqrt(mom.x() * mom.x() + mom.y() * mom.y()) / mom.z()));
1382   hitParams.phiSensX = (xMomentum == zMomentum ? phiSensX : -phiSensX);
1383   hitParams.phiSensY = (yMomentum == zMomentum ? phiSensY : -phiSensY);
1384 
1385   if (!hit.isValid()) {
1386     hitParams.hitState = TrackStruct::invalid;
1387     return hitParams;
1388   }
1389 
1390   // Get local positions and errors of hit and track
1391   const LocalPoint& lPHit = hit.localPosition();
1392   const LocalPoint& lPTrk = tsos.localPosition();
1393 
1394   // use APE also for the hit error, while APE is automatically included in tsos error
1395   //
1396   //  no need to add  APE to hitError anymore by Ajay 27 Oct 2014
1397   const LocalError& errHitApe = hit.localPositionError();  // now sum of CPE+APE as said by MARCO?
1398   LocalError errorWithoutAPE;
1399 
1400   bool Pixel(false);
1401   bool Strip(false);
1402 
1403   if (m_tkTreeVar_[rawId].subdetId == PixelSubdetector::PixelBarrel ||
1404       m_tkTreeVar_[rawId].subdetId == PixelSubdetector::PixelEndcap) {
1405     Pixel = true;
1406   } else if (m_tkTreeVar_[rawId].subdetId == StripSubdetector::TIB ||
1407              m_tkTreeVar_[rawId].subdetId == StripSubdetector::TOB ||
1408              m_tkTreeVar_[rawId].subdetId == StripSubdetector::TID ||
1409              m_tkTreeVar_[rawId].subdetId == StripSubdetector::TEC) {
1410     Strip = true;
1411   } else {
1412     edm::LogWarning("FillHitVariables") << "cant identify wether hit is from pixel or strip";
1413     hitParams.hitState = TrackStruct::invalid;
1414     return hitParams;
1415   }
1416 
1417   if (!hit.detUnit()) {
1418     hitParams.hitState = TrackStruct::invalid;
1419     return hitParams;
1420   }  // is it a single physical module?
1421   const GeomDetUnit& detUnit = *hit.detUnit();
1422 
1423   if (Pixel) {
1424     if (!dynamic_cast<const PixelTopology*>(&detUnit.type().topology())) {
1425       hitParams.hitState = TrackStruct::invalid;
1426       return hitParams;
1427     }
1428     const PixelGeomDetUnit* pixelDet = (const PixelGeomDetUnit*)(&detUnit);
1429     const LocalError& lape = pixelDet->localAlignmentError();
1430     if (lape.valid()) {
1431       errorWithoutAPE = LocalError(errHitApe.xx() - lape.xx(), errHitApe.xy() - lape.xy(), errHitApe.yy() - lape.yy());
1432     }
1433   }
1434   if (Strip) {
1435     if (!dynamic_cast<const StripTopology*>(&detUnit.type().topology())) {
1436       hitParams.hitState = TrackStruct::invalid;
1437       return hitParams;
1438     }
1439     const StripGeomDetUnit* stripDet = (const StripGeomDetUnit*)(&detUnit);
1440     const LocalError& lape = stripDet->localAlignmentError();
1441     if (lape.valid()) {
1442       errorWithoutAPE = LocalError(errHitApe.xx() - lape.xx(), errHitApe.xy() - lape.xy(), errHitApe.yy() - lape.yy());
1443     }
1444   }
1445 
1446   const LocalError& errHitWoApe = errorWithoutAPE;
1447   const LocalError& errTrk = tsos.localError().positionError();
1448 
1449   const StatePositionAndError2 positionAndError2Hit = this->positionAndError2(lPHit, errHitApe, hit);
1450   const StatePositionAndError2 positionAndError2HitWoApe = this->positionAndError2(lPHit, errHitWoApe, hit);
1451   edm::LogInfo("CalculateAPE") << "errHitWoApe  " << errHitWoApe << "errHitApe  " << errHitApe;
1452 
1453   const StatePositionAndError2 positionAndError2Trk = this->positionAndError2(lPTrk, errTrk, hit);
1454 
1455   const TrackStruct::HitState& stateHit(positionAndError2Hit.first);
1456   const TrackStruct::HitState& stateHitWoApe(positionAndError2HitWoApe.first);
1457   const TrackStruct::HitState& stateTrk(positionAndError2Trk.first);
1458 
1459   if (stateHit == TrackStruct::invalid || stateHitWoApe == TrackStruct::invalid || stateTrk == TrackStruct::invalid) {
1460     hitParams.hitState = TrackStruct::invalid;
1461     return hitParams;
1462   } else if (stateHit == TrackStruct::negativeError || stateHitWoApe == TrackStruct::negativeError ||
1463              stateTrk == TrackStruct::negativeError) {
1464     ++counter1;
1465     // Do not print error message by default
1466     //std::stringstream ss_error;
1467     //ss_error<<"Upper values belong to: ";
1468     //if(stateHit==TrackStruct::negativeError)ss_error<<"Hit without APE, ";
1469     //if(stateHitWoApe==TrackStruct::negativeError)ss_error<<"Hit with APE, ";
1470     //if(stateTrk==TrackStruct::negativeError)ss_error<<"Track,";
1471     //edm::LogError("Negative error Value")<<"@SUB=ApeEstimator::fillHitVariables"<<ss_error.str();
1472     hitParams.hitState = TrackStruct::negativeError;
1473     return hitParams;
1474   }
1475 
1476   // Calculate residuals
1477 
1478   const float xHit = positionAndError2Hit.second.posX;
1479   const float xTrk = positionAndError2Trk.second.posX;
1480   const float yHit = positionAndError2Hit.second.posY;
1481   const float yTrk = positionAndError2Trk.second.posY;
1482 
1483   const float errXHit2(positionAndError2Hit.second.errX2);
1484   const float errXHitWoApe2(positionAndError2HitWoApe.second.errX2);
1485   const float errXTrk2(positionAndError2Trk.second.errX2);
1486   const float errYHit2(positionAndError2Hit.second.errY2);
1487   const float errYHitWoApe2(positionAndError2HitWoApe.second.errY2);
1488   const float errYTrk2(positionAndError2Trk.second.errY2);
1489 
1490   const float errXHit = std::sqrt(positionAndError2Hit.second.errX2);
1491   const float errXHitWoApe = std::sqrt(positionAndError2HitWoApe.second.errX2);
1492   const float errXTrk = std::sqrt(positionAndError2Trk.second.errX2);
1493   const float errYHit = std::sqrt(positionAndError2Hit.second.errY2);
1494   const float errYHitWoApe = std::sqrt(positionAndError2HitWoApe.second.errY2);
1495   const float errYTrk = std::sqrt(positionAndError2Trk.second.errY2);
1496 
1497   const float resX = xTrk - xHit;
1498   const float resY = yTrk - yHit;
1499 
1500   const float errX = std::sqrt(errXHit2 + errXTrk2);
1501   const float errXWoApe2 = errXHitWoApe2 + errXTrk2;
1502   const float errXWoApe = std::sqrt(errXWoApe2);
1503   const float errY = std::sqrt(errYHit2 + errYTrk2);
1504   const float errYWoApe2 = errYHitWoApe2 + errYTrk2;
1505   const float errYWoApe = std::sqrt(errYWoApe2);
1506 
1507   const float norResX = resX / errX;
1508   const float norResY = resY / errY;
1509 
1510   // Take global orientation into account for residuals (sign is not important for errors)
1511 
1512   float resXprime(999.F), resYprime(999.F), norResXprime(999.F), norResYprime(999.F);
1513   if (m_tkTreeVar_[rawId].uDirection == 1) {
1514     resXprime = resX;
1515     norResXprime = norResX;
1516   } else if (m_tkTreeVar_[rawId].uDirection == -1) {
1517     resXprime = -resX;
1518     norResXprime = -norResX;
1519   } else {
1520     edm::LogError("FillHitVariables") << "Incorrect value of uDirection, which gives global module orientation";
1521     hitParams.hitState = TrackStruct::invalid;
1522     return hitParams;
1523   }
1524   if (m_tkTreeVar_[rawId].vDirection == 1) {
1525     resYprime = resY;
1526     norResYprime = norResY;
1527   } else if (m_tkTreeVar_[rawId].vDirection == -1) {
1528     resYprime = -resY;
1529     norResYprime = -norResY;
1530   } else {
1531     edm::LogError("FillHitVariables") << "Incorrect value of vDirection, which gives global module orientation";
1532     hitParams.hitState = TrackStruct::invalid;
1533     return hitParams;
1534   }
1535 
1536   hitParams.xHit = xHit;
1537   hitParams.xTrk = xTrk;
1538 
1539   hitParams.errXHit = errXHit;
1540   hitParams.errXHitWoApe = errXHitWoApe;
1541   hitParams.errXTrk = errXTrk;
1542 
1543   hitParams.errX2 = errX * errX;
1544   hitParams.errX = errX;
1545   hitParams.errXWoApe = errXWoApe;
1546 
1547   hitParams.resX = resXprime;
1548   hitParams.norResX = norResXprime;
1549 
1550   const float norResX2(norResXprime * norResXprime);
1551   hitParams.probX = TMath::Prob(norResX2, 1);
1552 
1553   hitParams.yHit = yHit;
1554   hitParams.yTrk = yTrk;
1555 
1556   hitParams.errYHit = errYHit;
1557   hitParams.errYHitWoApe = errYHitWoApe;
1558   hitParams.errYTrk = errYTrk;
1559 
1560   hitParams.errY2 = errY * errY;
1561   hitParams.errY = errY;
1562   hitParams.errYWoApe = errYWoApe;
1563 
1564   hitParams.resY = resYprime;
1565   hitParams.norResY = norResYprime;
1566 
1567   const float norResY2(norResYprime * norResYprime);
1568   hitParams.probY = TMath::Prob(norResY2, 1);
1569 
1570   // Cluster parameters
1571 
1572   if (m_tkTreeVar_[rawId].subdetId == PixelSubdetector::PixelBarrel ||
1573       m_tkTreeVar_[rawId].subdetId == PixelSubdetector::PixelEndcap) {
1574     const SiPixelRecHit& pixelHit = dynamic_cast<const SiPixelRecHit&>(recHit);
1575     const SiPixelCluster& pixelCluster = *pixelHit.cluster();
1576 
1577     hitParams.chargePixel = pixelCluster.charge();
1578     hitParams.widthX = pixelCluster.sizeX();
1579     hitParams.baryStripX = pixelCluster.x();
1580     hitParams.widthY = pixelCluster.sizeY();
1581     hitParams.baryStripY = pixelCluster.y();
1582 
1583     hitParams.clusterProbabilityXY = pixelHit.clusterProbability(0);
1584     hitParams.clusterProbabilityQ = pixelHit.clusterProbability(2);
1585     hitParams.clusterProbabilityXYQ = pixelHit.clusterProbability(1);
1586     hitParams.logClusterProbability = std::log10(hitParams.clusterProbabilityXY);
1587 
1588     hitParams.isOnEdge = pixelHit.isOnEdge();
1589     hitParams.hasBadPixels = pixelHit.hasBadPixels();
1590     hitParams.spansTwoRoc = pixelHit.spansTwoROCs();
1591     hitParams.qBin = pixelHit.qBin();
1592 
1593     hitParams.isPixelHit = true;
1594   } else if (m_tkTreeVar_[rawId].subdetId == StripSubdetector::TIB ||
1595              m_tkTreeVar_[rawId].subdetId == StripSubdetector::TOB ||
1596              m_tkTreeVar_[rawId].subdetId == StripSubdetector::TID ||
1597              m_tkTreeVar_[rawId].subdetId == StripSubdetector::TEC) {
1598     if (!(dynamic_cast<const SiStripRecHit2D*>(&recHit) || dynamic_cast<const SiStripRecHit1D*>(&recHit))) {
1599       edm::LogError("FillHitVariables")
1600           << "RecHit in Strip is 'Matched' or 'Projected', but here all should be monohits per module";
1601       hitParams.hitState = TrackStruct::invalid;
1602       return hitParams;
1603     }
1604     const SiStripCluster* clusterPtr(nullptr);
1605     if (m_tkTreeVar_[rawId].subdetId == StripSubdetector::TIB ||
1606         m_tkTreeVar_[rawId].subdetId == StripSubdetector::TOB) {
1607       if (dynamic_cast<const SiStripRecHit1D*>(&recHit)) {
1608         const SiStripRecHit1D& stripHit = dynamic_cast<const SiStripRecHit1D&>(recHit);
1609         clusterPtr = &(*stripHit.cluster());
1610       } else if (dynamic_cast<const SiStripRecHit2D*>(&recHit)) {
1611         edm::LogWarning("FillHitVariables") << "Data has TIB/TOB hits as SiStripRecHit2D and not 1D. Probably data is "
1612                                                "processed with CMSSW<34X. Nevertheless everything should work fine";
1613         const SiStripRecHit2D& stripHit = dynamic_cast<const SiStripRecHit2D&>(recHit);
1614         clusterPtr = &(*stripHit.cluster());
1615       }
1616     } else if (m_tkTreeVar_[rawId].subdetId == StripSubdetector::TID ||
1617                m_tkTreeVar_[rawId].subdetId == StripSubdetector::TEC) {
1618       const SiStripRecHit2D& stripHit = dynamic_cast<const SiStripRecHit2D&>(recHit);
1619       clusterPtr = &(*stripHit.cluster());
1620     }
1621     if (!clusterPtr) {
1622       edm::LogError("FillHitVariables") << "Pointer to cluster not valid!!! This should never happen...";
1623       hitParams.hitState = TrackStruct::invalid;
1624       return hitParams;
1625     }
1626     const SiStripCluster& stripCluster(*clusterPtr);
1627 
1628     siStripClusterInfo_.setCluster(stripCluster, rawId);
1629 
1630     auto stripChargeL = siStripClusterInfo_.stripCharges().begin();
1631     auto stripChargeR = siStripClusterInfo_.stripCharges().end() - 1;
1632     const std::pair<uint16_t, uint16_t> stripChargeLR = std::make_pair(*stripChargeL, *stripChargeR);
1633 
1634     hitParams.chargeStrip = siStripClusterInfo_.charge();
1635     hitParams.widthX = siStripClusterInfo_.width();
1636     hitParams.baryStripX = siStripClusterInfo_.baryStrip() + 1.;
1637     hitParams.isModuleUsable = siStripClusterInfo_.IsModuleUsable();
1638     hitParams.maxStrip = siStripClusterInfo_.maxStrip() + 1;
1639     hitParams.maxStripInv = m_tkTreeVar_[rawId].nStrips - hitParams.maxStrip + 1;
1640     hitParams.maxCharge = siStripClusterInfo_.maxCharge();
1641     hitParams.maxIndex = siStripClusterInfo_.maxIndex();
1642     hitParams.chargeOnEdges =
1643         static_cast<float>(stripChargeLR.first + stripChargeLR.second) / static_cast<float>(hitParams.chargeStrip);
1644     hitParams.chargeAsymmetry = static_cast<float>(stripChargeLR.first - stripChargeLR.second) /
1645                                 static_cast<float>(stripChargeLR.first + stripChargeLR.second);
1646     hitParams.chargeLRplus =
1647         static_cast<float>(siStripClusterInfo_.chargeLR().first + siStripClusterInfo_.chargeLR().second) /
1648         static_cast<float>(hitParams.chargeStrip);
1649     hitParams.chargeLRminus =
1650         static_cast<float>(siStripClusterInfo_.chargeLR().first - siStripClusterInfo_.chargeLR().second) /
1651         static_cast<float>(hitParams.chargeStrip);
1652     hitParams.sOverN = siStripClusterInfo_.signalOverNoise();
1653 
1654     // Calculate projection length corrected by drift
1655     if (!hit.detUnit()) {
1656       hitParams.hitState = TrackStruct::invalid;
1657       return hitParams;
1658     }  // is it a single physical module?
1659     const GeomDetUnit& detUnit = *hit.detUnit();
1660     if (!dynamic_cast<const StripTopology*>(&detUnit.type().topology())) {
1661       hitParams.hitState = TrackStruct::invalid;
1662       return hitParams;
1663     }
1664 
1665     const MagneticField* magField = &iSetup.getData(magFieldToken_);
1666     const SiStripLorentzAngle* lorentzAngle = &iSetup.getData(lorentzAngleToken_);
1667 
1668     const StripGeomDetUnit* stripDet = (const StripGeomDetUnit*)(&detUnit);
1669     LocalVector bField = (stripDet->surface()).toLocal(magField->inTesla(stripDet->surface().position()));
1670     float tanLorentzAnglePerTesla = lorentzAngle->getLorentzAngle(stripDet->geographicalId().rawId());
1671 
1672     float dirX = -tanLorentzAnglePerTesla * bField.y();
1673     float dirY = tanLorentzAnglePerTesla * bField.x();
1674     float dirZ = 1.;  // E field always in z direction
1675     LocalVector driftDirection(dirX, dirY, dirZ);
1676 
1677     const Bounds& bounds = stripDet->specificSurface().bounds();
1678     float maxLength = std::sqrt(std::pow(bounds.length(), 2) + std::pow(bounds.width(), 2));
1679     float thickness = bounds.thickness();
1680 
1681     const StripTopology& topol = dynamic_cast<const StripTopology&>(detUnit.type().topology());
1682     LocalVector momentumDir(tsos.localDirection());
1683     LocalPoint momentumPos(tsos.localPosition());
1684     LocalVector scaledMomentumDir(momentumDir);
1685     if (momentumDir.z() > 0.)
1686       scaledMomentumDir *= std::fabs(thickness / momentumDir.z());
1687     else if (momentumDir.z() < 0.)
1688       scaledMomentumDir *= -std::fabs(thickness / momentumDir.z());
1689     else
1690       scaledMomentumDir *= maxLength / momentumDir.mag();
1691 
1692     float projEdge1 = topol.measurementPosition(momentumPos - 0.5 * scaledMomentumDir).x();
1693     if (projEdge1 < 0.)
1694       projEdge1 = 0.;
1695     else if (projEdge1 > m_tkTreeVar_[rawId].nStrips)
1696       projEdge1 = m_tkTreeVar_[rawId].nStrips;
1697     float projEdge2 = topol.measurementPosition(momentumPos + 0.5 * scaledMomentumDir).x();
1698     if (projEdge2 < 0.)
1699       projEdge1 = 0.;
1700     else if (projEdge2 > m_tkTreeVar_[rawId].nStrips)
1701       projEdge1 = m_tkTreeVar_[rawId].nStrips;
1702 
1703     float coveredStrips = std::fabs(projEdge2 - projEdge1);
1704 
1705     hitParams.projWidth = coveredStrips;
1706 
1707   } else {
1708     edm::LogError("FillHitVariables") << "Incorrect subdetector ID, hit not associated to tracker";
1709     hitParams.hitState = TrackStruct::notInTracker;
1710     return hitParams;
1711   }
1712 
1713   if (!hitParams.isModuleUsable) {
1714     hitParams.hitState = TrackStruct::invalid;
1715     return hitParams;
1716   }
1717 
1718   if (hitParams.v_sector.empty()) {
1719     hitParams.hitState = TrackStruct::notAssignedToSectors;
1720     return hitParams;
1721   }
1722 
1723   return hitParams;
1724   //}
1725 }
1726 
1727 ApeEstimator::StatePositionAndError2 ApeEstimator::positionAndError2(const LocalPoint& localPoint,
1728                                                                      const LocalError& localError,
1729                                                                      const TransientTrackingRecHit& hit) {
1730   StatePositionAndError2 vPE2 = std::make_pair(TrackStruct::invalid, PositionAndError2());
1731 
1732   const DetId& detId(hit.geographicalId());
1733   const uint32_t& rawId(detId.rawId());
1734   const unsigned int& subdetId(m_tkTreeVar_[rawId].subdetId);
1735 
1736   if (localError.xx() < 0. || localError.yy() < 0.) {
1737     // Do not print error message by default
1738     //edm::LogError("Negative error Value")<<"@SUB=ApeEstimator::fillHitVariables"
1739     //                                     <<"One of the squared error methods gives negative result\n"
1740     //                                     <<"\tSubdetector\tlocalError.xx()\tlocalError.yy()\n"
1741     //                                     <<"\t"<<subdetId<<"\t\t"<<localError.xx()<<"\t"<<localError.yy();
1742     vPE2.first = TrackStruct::negativeError;
1743     return vPE2;
1744   }
1745 
1746   if (subdetId == PixelSubdetector::PixelBarrel || subdetId == PixelSubdetector::PixelEndcap ||
1747       subdetId == StripSubdetector::TIB || subdetId == StripSubdetector::TOB) {
1748     // Cartesian coordinates
1749     vPE2 = std::make_pair(TrackStruct::ok, this->rectangularPositionAndError2(localPoint, localError));
1750   } else if (subdetId == StripSubdetector::TID || subdetId == StripSubdetector::TEC) {
1751     // Local x in radial coordinates
1752     if (!hit.detUnit())
1753       return vPE2;  // is it a single physical module?
1754     const GeomDetUnit& detUnit = *hit.detUnit();
1755 
1756     if (!dynamic_cast<const RadialStripTopology*>(&detUnit.type().topology()))
1757       return vPE2;
1758     const RadialStripTopology& topol = dynamic_cast<const RadialStripTopology&>(detUnit.type().topology());
1759 
1760     MeasurementError measError = topol.measurementError(localPoint, localError);
1761     if (measError.uu() < 0. || measError.vv() < 0.) {
1762       // Do not print error message by default
1763       //edm::LogError("Negative error Value")<<"@SUB=ApeEstimator::fillHitVariables"
1764       //                                     <<"One of the squared error methods gives negative result\n"
1765       //                                     <<"\tmeasError.uu()\tmeasError.vv()\n"
1766       //                                     <<"\t"<<measError.uu()<<"\t"<<measError.vv()
1767       //                                     <<"\n\nOriginalValues:\n"
1768       //                                     <<localPoint.x()<<" "<<localPoint.y()<<"\n"
1769       //                                     <<localError.xx()<<" "<<localError.yy()<<"\n"
1770       //                                     <<"Subdet: "<<subdetId;
1771       vPE2.first = TrackStruct::negativeError;
1772       return vPE2;
1773     }
1774     vPE2 = std::make_pair(TrackStruct::ok, this->radialPositionAndError2(localPoint, localError, topol));
1775   } else {
1776     edm::LogError("FillHitVariables") << "Incorrect subdetector ID, hit not associated to tracker";
1777   }
1778 
1779   return vPE2;
1780 }
1781 
1782 ApeEstimator::PositionAndError2 ApeEstimator::rectangularPositionAndError2(const LocalPoint& lP, const LocalError& lE) {
1783   const float x(lP.x());
1784   const float y(lP.y());
1785   const float errX2(lE.xx());
1786   const float errY2(lE.yy());
1787 
1788   return PositionAndError2(x, y, errX2, errY2);
1789 }
1790 
1791 ApeEstimator::PositionAndError2 ApeEstimator::radialPositionAndError2(const LocalPoint& lP,
1792                                                                       const LocalError& lE,
1793                                                                       const RadialStripTopology& topol) {
1794   MeasurementPoint measPos = topol.measurementPosition(lP);
1795   MeasurementError measErr = topol.measurementError(lP, lE);
1796 
1797   const float r_0 = topol.originToIntersection();
1798   const float stripLength = topol.localStripLength(lP);
1799   const float phi = topol.stripAngle(measPos.x());
1800 
1801   float x(-999.F);
1802   float y(-999.F);
1803   float errX2(-999.F);
1804   float errY2(-999.F);
1805 
1806   x = phi * r_0;
1807   // Radial y (not symmetric around 0; radial distance with minimum at middle strip at lower edge [0, yMax])
1808   const float l_0 = r_0 - topol.detHeight() / 2;
1809   const float cosPhi(std::cos(phi));
1810   y = measPos.y() * stripLength - 0.5 * stripLength + l_0 * (1. / cosPhi - 1.);
1811 
1812   const float angularWidth2(topol.angularWidth() * topol.angularWidth());
1813   const float errPhi2(measErr.uu() * angularWidth2);
1814 
1815   errX2 = errPhi2 * r_0 * r_0;
1816   // Radial y (not symmetric around 0, real radial distance from intersection point)
1817   const float cosPhi4(std::pow(cosPhi, 4)), sinPhi2(std::sin(phi) * std::sin(phi));
1818   const float helpSummand = l_0 * l_0 * (sinPhi2 / cosPhi4 * errPhi2);
1819   errY2 = measErr.vv() * stripLength * stripLength + helpSummand;
1820 
1821   return PositionAndError2(x, y, errX2, errY2);
1822 }
1823 
1824 // -----------------------------------------------------------------------------------------------------------
1825 
1826 void ApeEstimator::hitSelection() {
1827   this->setHitSelectionMapUInt("width");
1828   this->setHitSelectionMap("widthProj");
1829   this->setHitSelectionMap("widthDiff");
1830   this->setHitSelectionMap("charge");
1831   this->setHitSelectionMapUInt("edgeStrips");
1832   this->setHitSelectionMap("maxCharge");
1833   this->setHitSelectionMapUInt("maxIndex");
1834   this->setHitSelectionMap("chargeOnEdges");
1835   this->setHitSelectionMap("chargeAsymmetry");
1836   this->setHitSelectionMap("chargeLRplus");
1837   this->setHitSelectionMap("chargeLRminus");
1838   this->setHitSelectionMap("sOverN");
1839 
1840   this->setHitSelectionMap("chargePixel");
1841   this->setHitSelectionMapUInt("widthX");
1842   this->setHitSelectionMapUInt("widthY");
1843 
1844   this->setHitSelectionMap("baryStripX");
1845   this->setHitSelectionMap("baryStripY");
1846   this->setHitSelectionMap("clusterProbabilityXY");
1847   this->setHitSelectionMap("clusterProbabilityQ");
1848   this->setHitSelectionMap("clusterProbabilityXYQ");
1849   this->setHitSelectionMap("logClusterProbability");
1850   this->setHitSelectionMapUInt("isOnEdge");
1851   this->setHitSelectionMapUInt("hasBadPixels");
1852   this->setHitSelectionMapUInt("spansTwoRoc");
1853   this->setHitSelectionMapUInt("qBin");
1854 
1855   this->setHitSelectionMap("phiSens");
1856   this->setHitSelectionMap("phiSensX");
1857   this->setHitSelectionMap("phiSensY");
1858   this->setHitSelectionMap("resX");
1859   this->setHitSelectionMap("norResX");
1860   this->setHitSelectionMap("probX");
1861   this->setHitSelectionMap("errXHit");
1862   this->setHitSelectionMap("errXTrk");
1863   this->setHitSelectionMap("errX");
1864   this->setHitSelectionMap("errX2");
1865 
1866   this->setHitSelectionMap("resY");
1867   this->setHitSelectionMap("norResY");
1868   this->setHitSelectionMap("probY");
1869   this->setHitSelectionMap("errYHit");
1870   this->setHitSelectionMap("errYTrk");
1871   this->setHitSelectionMap("errY");
1872   this->setHitSelectionMap("errY2");
1873 
1874   edm::LogInfo("HitSelector") << "applying hit cuts ...";
1875   bool emptyMap(true);
1876   for (auto& i_hitSelection : m_hitSelection_) {
1877     if (!i_hitSelection.second.empty()) {
1878       int entry(1);
1879       double intervalBegin(999.);
1880       for (std::vector<double>::iterator i_hitInterval = i_hitSelection.second.begin();
1881            i_hitInterval != i_hitSelection.second.end();
1882            ++entry) {
1883         if (entry % 2 == 1) {
1884           intervalBegin = *i_hitInterval;
1885           ++i_hitInterval;
1886         } else {
1887           if (intervalBegin > *i_hitInterval) {
1888             edm::LogError("HitSelector") << "INVALID Interval selected for  " << i_hitSelection.first << ":\t"
1889                                          << intervalBegin << " > " << (*i_hitInterval) << "\n ... delete Selection for "
1890                                          << i_hitSelection.first;
1891             i_hitSelection.second.clear();
1892             i_hitInterval = i_hitSelection.second.begin();  //emptyMap = true; i_hitSelection = m_hitSelection_.begin();
1893           } else {
1894             edm::LogInfo("HitSelector") << "Interval selected for  " << i_hitSelection.first << ":\t" << intervalBegin
1895                                         << ", " << (*i_hitInterval);
1896             ++i_hitInterval;
1897           }
1898         }
1899       }
1900       if (!i_hitSelection.second.empty())
1901         emptyMap = false;
1902     }
1903   }
1904 
1905   bool emptyMapUInt(true);
1906   for (auto& i_hitSelection : m_hitSelectionUInt_) {
1907     if (!i_hitSelection.second.empty()) {
1908       int entry(1);
1909       unsigned int intervalBegin(999);
1910       for (std::vector<unsigned int>::iterator i_hitInterval = i_hitSelection.second.begin();
1911            i_hitInterval != i_hitSelection.second.end();
1912            ++entry) {
1913         if (entry % 2 == 1) {
1914           intervalBegin = *i_hitInterval;
1915           ++i_hitInterval;
1916         } else {
1917           if (intervalBegin > *i_hitInterval) {
1918             edm::LogError("HitSelector") << "INVALID Interval selected for  " << i_hitSelection.first << ":\t"
1919                                          << intervalBegin << " > " << (*i_hitInterval) << "\n ... delete Selection for "
1920                                          << i_hitSelection.first;
1921             i_hitSelection.second.clear();
1922             i_hitInterval = i_hitSelection.second.begin();  //emptyMap = true; i_hitSelection = m_hitSelection_.begin();
1923           } else {
1924             edm::LogInfo("HitSelector") << "Interval selected for  " << i_hitSelection.first << ":\t" << intervalBegin
1925                                         << ", " << (*i_hitInterval);
1926             ++i_hitInterval;
1927           }
1928         }
1929       }
1930       if (!i_hitSelection.second.empty())
1931         emptyMapUInt = false;
1932     }
1933   }
1934 
1935   if (emptyMap && emptyMapUInt) {
1936     m_hitSelection_.clear();
1937     m_hitSelectionUInt_.clear();
1938     edm::LogInfo("HitSelector") << "NO hit cuts applied";
1939   }
1940   return;
1941 }
1942 
1943 void ApeEstimator::setHitSelectionMap(const std::string& cutVariable) {
1944   edm::ParameterSet parSet(parameterSet_.getParameter<edm::ParameterSet>("HitSelector"));
1945   std::vector<double> v_cutVariable(parSet.getParameter<std::vector<double> >(cutVariable));
1946   if (v_cutVariable.size() % 2 == 1) {
1947     edm::LogError("HitSelector") << "Invalid Hit Selection for " << cutVariable
1948                                  << ": need even number of arguments (intervals)"
1949                                  << "\n ... delete Selection for " << cutVariable;
1950     v_cutVariable.clear();
1951     m_hitSelection_[cutVariable] = v_cutVariable;
1952     return;
1953   }
1954   m_hitSelection_[cutVariable] = v_cutVariable;
1955   return;
1956 }
1957 
1958 void ApeEstimator::setHitSelectionMapUInt(const std::string& cutVariable) {
1959   edm::ParameterSet parSet(parameterSet_.getParameter<edm::ParameterSet>("HitSelector"));
1960   std::vector<unsigned int> v_cutVariable(parSet.getParameter<std::vector<unsigned int> >(cutVariable));
1961   if (v_cutVariable.size() % 2 == 1) {
1962     edm::LogError("HitSelector") << "Invalid Hit Selection for " << cutVariable
1963                                  << ": need even number of arguments (intervals)"
1964                                  << "\n ... delete Selection for " << cutVariable;
1965     v_cutVariable.clear();
1966     m_hitSelectionUInt_[cutVariable] = v_cutVariable;
1967     return;
1968   }
1969   m_hitSelectionUInt_[cutVariable] = v_cutVariable;
1970   return;
1971 }
1972 
1973 // -----------------------------------------------------------------------------------------------------------
1974 
1975 bool ApeEstimator::hitSelected(TrackStruct::HitParameterStruct& hitParams) const {
1976   if (hitParams.hitState == TrackStruct::notInTracker)
1977     return false;
1978   if (hitParams.hitState == TrackStruct::invalid || hitParams.hitState == TrackStruct::negativeError)
1979     return false;
1980 
1981   bool isGoodHit(true);
1982   bool isGoodHitX(true);
1983   bool isGoodHitY(true);
1984 
1985   for (auto& i_hitSelection : m_hitSelection_) {
1986     const std::string& hitSelection(i_hitSelection.first);
1987     const std::vector<double>& v_hitSelection(i_hitSelection.second);
1988     if (v_hitSelection.empty())
1989       continue;
1990 
1991     // For pixel and strip sectors in common
1992     if (hitSelection == "phiSens") {
1993       if (!this->inDoubleInterval(v_hitSelection, hitParams.phiSens))
1994         isGoodHit = false;
1995     } else if (hitSelection == "phiSensX") {
1996       if (!this->inDoubleInterval(v_hitSelection, hitParams.phiSensX))
1997         isGoodHit = false;
1998     } else if (hitSelection == "phiSensY") {
1999       if (!this->inDoubleInterval(v_hitSelection, hitParams.phiSensY))
2000         isGoodHit = false;
2001     }
2002 
2003     else if (hitSelection == "resX") {
2004       if (!this->inDoubleInterval(v_hitSelection, hitParams.resX))
2005         isGoodHitX = false;
2006     } else if (hitSelection == "norResX") {
2007       if (!this->inDoubleInterval(v_hitSelection, hitParams.norResX))
2008         isGoodHitX = false;
2009     } else if (hitSelection == "probX") {
2010       if (!this->inDoubleInterval(v_hitSelection, hitParams.probX))
2011         isGoodHitX = false;
2012     } else if (hitSelection == "errXHit") {
2013       if (!this->inDoubleInterval(v_hitSelection, hitParams.errXHit))
2014         isGoodHitX = false;
2015     } else if (hitSelection == "errXTrk") {
2016       if (!this->inDoubleInterval(v_hitSelection, hitParams.errXTrk))
2017         isGoodHitX = false;
2018     } else if (hitSelection == "errX") {
2019       if (!this->inDoubleInterval(v_hitSelection, hitParams.errX))
2020         isGoodHitX = false;
2021     } else if (hitSelection == "errX2") {
2022       if (!this->inDoubleInterval(v_hitSelection, hitParams.errX2))
2023         isGoodHitX = false;
2024     }
2025 
2026     // For pixel only
2027     if (hitParams.isPixelHit) {
2028       if (hitSelection == "chargePixel") {
2029         if (!this->inDoubleInterval(v_hitSelection, hitParams.chargePixel))
2030           isGoodHit = false;
2031       } else if (hitSelection == "clusterProbabilityXY") {
2032         if (!this->inDoubleInterval(v_hitSelection, hitParams.clusterProbabilityXY))
2033           isGoodHit = false;
2034       } else if (hitSelection == "clusterProbabilityQ") {
2035         if (!this->inDoubleInterval(v_hitSelection, hitParams.clusterProbabilityQ))
2036           isGoodHit = false;
2037       } else if (hitSelection == "clusterProbabilityXYQ") {
2038         if (!this->inDoubleInterval(v_hitSelection, hitParams.clusterProbabilityXYQ))
2039           isGoodHit = false;
2040       } else if (hitSelection == "logClusterProbability") {
2041         if (!this->inDoubleInterval(v_hitSelection, hitParams.logClusterProbability))
2042           isGoodHit = false;
2043       }
2044 
2045       else if (hitSelection == "baryStripX") {
2046         if (!this->inDoubleInterval(v_hitSelection, hitParams.baryStripX))
2047           isGoodHitX = false;
2048       } else if (hitSelection == "baryStripY") {
2049         if (!this->inDoubleInterval(v_hitSelection, hitParams.baryStripY))
2050           isGoodHitY = false;
2051       }
2052 
2053       else if (hitSelection == "resY") {
2054         if (!this->inDoubleInterval(v_hitSelection, hitParams.resY))
2055           isGoodHitY = false;
2056       } else if (hitSelection == "norResY") {
2057         if (!this->inDoubleInterval(v_hitSelection, hitParams.norResY))
2058           isGoodHitY = false;
2059       } else if (hitSelection == "probY") {
2060         if (!this->inDoubleInterval(v_hitSelection, hitParams.probY))
2061           isGoodHitY = false;
2062       } else if (hitSelection == "errYHit") {
2063         if (!this->inDoubleInterval(v_hitSelection, hitParams.errYHit))
2064           isGoodHitY = false;
2065       } else if (hitSelection == "errYTrk") {
2066         if (!this->inDoubleInterval(v_hitSelection, hitParams.errYTrk))
2067           isGoodHitY = false;
2068       } else if (hitSelection == "errY") {
2069         if (!this->inDoubleInterval(v_hitSelection, hitParams.errY))
2070           isGoodHitY = false;
2071       } else if (hitSelection == "errY2") {
2072         if (!this->inDoubleInterval(v_hitSelection, hitParams.errY2))
2073           isGoodHitY = false;
2074       }
2075     } else {  // For strip only
2076       if (hitSelection == "widthProj") {
2077         if (!this->inDoubleInterval(v_hitSelection, hitParams.projWidth))
2078           isGoodHit = false;
2079       } else if (hitSelection == "widthDiff") {
2080         if (!this->inDoubleInterval(v_hitSelection, hitParams.projWidth - static_cast<float>(hitParams.widthX)))
2081           isGoodHit = false;
2082       } else if (hitSelection == "charge") {
2083         if (!this->inDoubleInterval(v_hitSelection, hitParams.chargeStrip))
2084           isGoodHit = false;
2085       } else if (hitSelection == "maxCharge") {
2086         if (!this->inDoubleInterval(v_hitSelection, hitParams.maxCharge))
2087           isGoodHit = false;
2088       } else if (hitSelection == "chargeOnEdges") {
2089         if (!this->inDoubleInterval(v_hitSelection, hitParams.chargeOnEdges))
2090           isGoodHit = false;
2091       } else if (hitSelection == "chargeAsymmetry") {
2092         if (!this->inDoubleInterval(v_hitSelection, hitParams.chargeAsymmetry))
2093           isGoodHit = false;
2094       } else if (hitSelection == "chargeLRplus") {
2095         if (!this->inDoubleInterval(v_hitSelection, hitParams.chargeLRplus))
2096           isGoodHit = false;
2097       } else if (hitSelection == "chargeLRminus") {
2098         if (!this->inDoubleInterval(v_hitSelection, hitParams.chargeLRminus))
2099           isGoodHit = false;
2100       } else if (hitSelection == "sOverN") {
2101         if (!this->inDoubleInterval(v_hitSelection, hitParams.sOverN))
2102           isGoodHit = false;
2103       }
2104     }
2105   }
2106 
2107   for (auto& i_hitSelection : m_hitSelectionUInt_) {
2108     const std::string& hitSelection(i_hitSelection.first);
2109     const std::vector<unsigned int>& v_hitSelection(i_hitSelection.second);
2110     if (v_hitSelection.empty())
2111       continue;
2112 
2113     // For pixel and strip sectors in common
2114 
2115     // For pixel only
2116     if (hitParams.isPixelHit) {
2117       if (hitSelection == "isOnEdge") {
2118         if (!this->inUintInterval(v_hitSelection, hitParams.isOnEdge))
2119           isGoodHit = false;
2120       } else if (hitSelection == "hasBadPixels") {
2121         if (!this->inUintInterval(v_hitSelection, hitParams.hasBadPixels))
2122           isGoodHit = false;
2123       } else if (hitSelection == "spansTwoRoc") {
2124         if (!this->inUintInterval(v_hitSelection, hitParams.spansTwoRoc))
2125           isGoodHit = false;
2126       } else if (hitSelection == "qBin") {
2127         if (!this->inUintInterval(v_hitSelection, hitParams.qBin))
2128           isGoodHit = false;
2129       }
2130 
2131       else if (hitSelection == "widthX") {
2132         if (!this->inUintInterval(v_hitSelection, hitParams.widthX))
2133           isGoodHitX = false;
2134       } else if (hitSelection == "widthY") {
2135         if (!this->inUintInterval(v_hitSelection, hitParams.widthY))
2136           isGoodHitY = false;
2137       }
2138     } else {  // For strip only
2139       if (hitSelection == "width") {
2140         if (!this->inUintInterval(v_hitSelection, hitParams.widthX))
2141           isGoodHit = false;
2142       } else if (hitSelection == "edgeStrips") {
2143         if (!this->inUintInterval(v_hitSelection, hitParams.maxStrip, hitParams.maxStripInv))
2144           isGoodHit = false;
2145       } else if (hitSelection == "maxIndex") {
2146         if (!this->inUintInterval(v_hitSelection, hitParams.maxIndex))
2147           isGoodHit = false;
2148       }
2149     }
2150   }
2151 
2152   if (hitParams.isPixelHit) {
2153     hitParams.goodXMeasurement = isGoodHit && isGoodHitX;
2154     hitParams.goodYMeasurement = isGoodHit && isGoodHitY;
2155   } else {
2156     hitParams.goodXMeasurement = isGoodHit && isGoodHitX;
2157     hitParams.goodYMeasurement = false;
2158   }
2159 
2160   if (!hitParams.goodXMeasurement && !hitParams.goodYMeasurement)
2161     return false;
2162   else
2163     return true;
2164 }
2165 
2166 bool ApeEstimator::inDoubleInterval(const std::vector<double>& v_hitSelection, const float variable) const {
2167   int entry(0);
2168   double intervalBegin(999.);
2169   bool isSelected(false);
2170   for (auto const& i_hitInterval : v_hitSelection) {
2171     ++entry;
2172     if (entry % 2 == 1)
2173       intervalBegin = i_hitInterval;
2174     else if (variable >= intervalBegin && variable < i_hitInterval)
2175       isSelected = true;
2176   }
2177   return isSelected;
2178 }
2179 
2180 bool ApeEstimator::inUintInterval(const std::vector<unsigned int>& v_hitSelection,
2181                                   const unsigned int variable,
2182                                   const unsigned int variable2) const {
2183   int entry(0);
2184   unsigned int intervalBegin(999);
2185   bool isSelected(false);
2186   for (auto i_hitInterval : v_hitSelection) {
2187     ++entry;
2188     if (entry % 2 == 1)
2189       intervalBegin = i_hitInterval;
2190     else if (variable >= intervalBegin && variable <= i_hitInterval) {
2191       if (variable2 == 999 || (variable2 >= intervalBegin && variable2 <= i_hitInterval))
2192         isSelected = true;
2193     }
2194   }
2195   return isSelected;
2196 }
2197 
2198 // -----------------------------------------------------------------------------------------------------------
2199 
2200 void ApeEstimator::fillHistsForAnalyzerMode(const TrackStruct& trackStruct) {
2201   unsigned int goodHitsPerTrack(trackStruct.v_hitParams.size());
2202   tkDetector_.HitsGood->Fill(goodHitsPerTrack);
2203   tkDetector_.HitsGoodVsHitsValid->Fill(trackStruct.trkParams.hitsValid, goodHitsPerTrack);
2204   tkDetector_.PHitsGoodVsHitsValid->Fill(trackStruct.trkParams.hitsValid, goodHitsPerTrack);
2205 
2206   if (parameterSet_.getParameter<bool>("applyTrackCuts")) {
2207     // which tracks to take? need min. nr. of selected hits?
2208     if (goodHitsPerTrack < minGoodHitsPerTrack_)
2209       return;
2210   }
2211 
2212   tkDetector_.HitsSize->Fill(trackStruct.trkParams.hitsSize);
2213   tkDetector_.HitsValid->Fill(trackStruct.trkParams.hitsValid);
2214   tkDetector_.HitsInvalid->Fill(trackStruct.trkParams.hitsInvalid);
2215   tkDetector_.Hits2D->Fill(trackStruct.trkParams.hits2D);
2216   tkDetector_.LayersMissed->Fill(trackStruct.trkParams.layersMissed);
2217   tkDetector_.HitsPixel->Fill(trackStruct.trkParams.hitsPixel);
2218   tkDetector_.HitsStrip->Fill(trackStruct.trkParams.hitsStrip);
2219   tkDetector_.Charge->Fill(trackStruct.trkParams.charge);
2220   tkDetector_.Chi2->Fill(trackStruct.trkParams.chi2);
2221   tkDetector_.Ndof->Fill(trackStruct.trkParams.ndof);
2222   tkDetector_.NorChi2->Fill(trackStruct.trkParams.norChi2);
2223   tkDetector_.Prob->Fill(trackStruct.trkParams.prob);
2224   tkDetector_.Eta->Fill(trackStruct.trkParams.eta);
2225   tkDetector_.EtaErr->Fill(trackStruct.trkParams.etaErr);
2226   tkDetector_.EtaSig->Fill(trackStruct.trkParams.eta / trackStruct.trkParams.etaErr);
2227   tkDetector_.Theta->Fill(trackStruct.trkParams.theta * 180. / M_PI);
2228   tkDetector_.Phi->Fill(trackStruct.trkParams.phi * 180. / M_PI);
2229   tkDetector_.PhiErr->Fill(trackStruct.trkParams.phiErr * 180. / M_PI);
2230   tkDetector_.PhiSig->Fill(trackStruct.trkParams.phi / trackStruct.trkParams.phiErr);
2231   tkDetector_.D0Beamspot->Fill(trackStruct.trkParams.d0Beamspot);
2232   tkDetector_.D0BeamspotErr->Fill(trackStruct.trkParams.d0BeamspotErr);
2233   tkDetector_.D0BeamspotSig->Fill(trackStruct.trkParams.d0Beamspot / trackStruct.trkParams.d0BeamspotErr);
2234   tkDetector_.Dz->Fill(trackStruct.trkParams.dz);
2235   tkDetector_.DzErr->Fill(trackStruct.trkParams.dzErr);
2236   tkDetector_.DzSig->Fill(trackStruct.trkParams.dz / trackStruct.trkParams.dzErr);
2237   tkDetector_.P->Fill(trackStruct.trkParams.p);
2238   tkDetector_.Pt->Fill(trackStruct.trkParams.pt);
2239   tkDetector_.PtErr->Fill(trackStruct.trkParams.ptErr);
2240   tkDetector_.PtSig->Fill(trackStruct.trkParams.pt / trackStruct.trkParams.ptErr);
2241   tkDetector_.MeanAngle->Fill(trackStruct.trkParams.meanPhiSensToNorm * 180. / M_PI);
2242 
2243   tkDetector_.MeanAngleVsHits->Fill(trackStruct.trkParams.hitsSize,
2244                                     trackStruct.trkParams.meanPhiSensToNorm * 180. / M_PI);
2245   tkDetector_.HitsPixelVsEta->Fill(trackStruct.trkParams.eta, trackStruct.trkParams.hitsPixel);
2246   tkDetector_.HitsPixelVsTheta->Fill(trackStruct.trkParams.theta * 180. / M_PI, trackStruct.trkParams.hitsPixel);
2247   tkDetector_.HitsStripVsEta->Fill(trackStruct.trkParams.eta, trackStruct.trkParams.hitsStrip);
2248   tkDetector_.HitsStripVsTheta->Fill(trackStruct.trkParams.theta * 180. / M_PI, trackStruct.trkParams.hitsStrip);
2249   tkDetector_.PtVsEta->Fill(trackStruct.trkParams.eta, trackStruct.trkParams.pt);
2250   tkDetector_.PtVsTheta->Fill(trackStruct.trkParams.theta * 180. / M_PI, trackStruct.trkParams.pt);
2251 
2252   tkDetector_.PMeanAngleVsHits->Fill(trackStruct.trkParams.hitsSize,
2253                                      trackStruct.trkParams.meanPhiSensToNorm * 180. / M_PI);
2254   tkDetector_.PHitsPixelVsEta->Fill(trackStruct.trkParams.eta, trackStruct.trkParams.hitsPixel);
2255   tkDetector_.PHitsPixelVsTheta->Fill(trackStruct.trkParams.theta * 180. / M_PI, trackStruct.trkParams.hitsPixel);
2256   tkDetector_.PHitsStripVsEta->Fill(trackStruct.trkParams.eta, trackStruct.trkParams.hitsStrip);
2257   tkDetector_.PHitsStripVsTheta->Fill(trackStruct.trkParams.theta * 180. / M_PI, trackStruct.trkParams.hitsStrip);
2258   tkDetector_.PPtVsEta->Fill(trackStruct.trkParams.eta, trackStruct.trkParams.pt);
2259   tkDetector_.PPtVsTheta->Fill(trackStruct.trkParams.theta * 180. / M_PI, trackStruct.trkParams.pt);
2260 
2261   for (auto& i_hit : trackStruct.v_hitParams) {
2262     const TrackStruct::HitParameterStruct& hit(i_hit);
2263     //Put here from earlier method
2264     if (hit.hitState == TrackStruct::notAssignedToSectors)
2265       continue;
2266 
2267     for (auto& i_sector : m_tkSector_) {
2268       bool moduleInSector(false);
2269       for (auto const& i_hitSector : hit.v_sector) {
2270         if (i_sector.first == i_hitSector) {
2271           moduleInSector = true;
2272           break;
2273         }
2274       }
2275       if (!moduleInSector)
2276         continue;
2277       TrackerSectorStruct& sector(i_sector.second);
2278 
2279       if (hit.goodXMeasurement) {
2280         std::map<std::string, TrackerSectorStruct::CorrelationHists>& m_corrHists(sector.m_correlationHistsX);
2281 
2282         // Cluster and Hit Parameters
2283         this->fillHitHistsXForAnalyzerMode(hit, sector);
2284 
2285         // Track Parameters
2286         m_corrHists["HitsValid"].fillCorrHistsX(hit, trackStruct.trkParams.hitsValid);
2287         m_corrHists["HitsGood"].fillCorrHistsX(hit, goodHitsPerTrack);
2288         m_corrHists["HitsInvalid"].fillCorrHistsX(hit, trackStruct.trkParams.hitsInvalid);
2289         m_corrHists["Hits2D"].fillCorrHistsX(hit, trackStruct.trkParams.hits2D);
2290         m_corrHists["LayersMissed"].fillCorrHistsX(hit, trackStruct.trkParams.layersMissed);
2291         m_corrHists["HitsPixel"].fillCorrHistsX(hit, trackStruct.trkParams.hitsPixel);
2292         m_corrHists["HitsStrip"].fillCorrHistsX(hit, trackStruct.trkParams.hitsStrip);
2293         m_corrHists["NorChi2"].fillCorrHistsX(hit, trackStruct.trkParams.norChi2);
2294         m_corrHists["Theta"].fillCorrHistsX(hit, trackStruct.trkParams.theta * 180. / M_PI);
2295         m_corrHists["Phi"].fillCorrHistsX(hit, trackStruct.trkParams.phi * 180. / M_PI);
2296         m_corrHists["D0Beamspot"].fillCorrHistsX(hit, trackStruct.trkParams.d0Beamspot);
2297         m_corrHists["Dz"].fillCorrHistsX(hit, trackStruct.trkParams.dz);
2298         m_corrHists["Pt"].fillCorrHistsX(hit, trackStruct.trkParams.pt);
2299         m_corrHists["P"].fillCorrHistsX(hit, trackStruct.trkParams.p);
2300         m_corrHists["InvP"].fillCorrHistsX(hit, 1. / trackStruct.trkParams.p);
2301         m_corrHists["MeanAngle"].fillCorrHistsX(hit, trackStruct.trkParams.meanPhiSensToNorm * 180. / M_PI);
2302       }
2303 
2304       if (hit.goodYMeasurement) {
2305         std::map<std::string, TrackerSectorStruct::CorrelationHists>& m_corrHists(sector.m_correlationHistsY);
2306 
2307         // Cluster and Hit Parameters
2308         this->fillHitHistsYForAnalyzerMode(hit, sector);
2309 
2310         // Track Parameters
2311         m_corrHists["HitsValid"].fillCorrHistsY(hit, trackStruct.trkParams.hitsValid);
2312         m_corrHists["HitsGood"].fillCorrHistsY(hit, goodHitsPerTrack);
2313         m_corrHists["HitsInvalid"].fillCorrHistsY(hit, trackStruct.trkParams.hitsInvalid);
2314         m_corrHists["Hits2D"].fillCorrHistsY(hit, trackStruct.trkParams.hits2D);
2315         m_corrHists["LayersMissed"].fillCorrHistsY(hit, trackStruct.trkParams.layersMissed);
2316         m_corrHists["HitsPixel"].fillCorrHistsY(hit, trackStruct.trkParams.hitsPixel);
2317         m_corrHists["HitsStrip"].fillCorrHistsY(hit, trackStruct.trkParams.hitsStrip);
2318         m_corrHists["NorChi2"].fillCorrHistsY(hit, trackStruct.trkParams.norChi2);
2319         m_corrHists["Theta"].fillCorrHistsY(hit, trackStruct.trkParams.theta * 180. / M_PI);
2320         m_corrHists["Phi"].fillCorrHistsY(hit, trackStruct.trkParams.phi * 180. / M_PI);
2321         m_corrHists["D0Beamspot"].fillCorrHistsY(hit, trackStruct.trkParams.d0Beamspot);
2322         m_corrHists["Dz"].fillCorrHistsY(hit, trackStruct.trkParams.dz);
2323         m_corrHists["Pt"].fillCorrHistsY(hit, trackStruct.trkParams.pt);
2324         m_corrHists["P"].fillCorrHistsY(hit, trackStruct.trkParams.p);
2325         m_corrHists["InvP"].fillCorrHistsY(hit, 1. / trackStruct.trkParams.p);
2326         m_corrHists["MeanAngle"].fillCorrHistsY(hit, trackStruct.trkParams.meanPhiSensToNorm * 180. / M_PI);
2327       }
2328 
2329       // Special Histograms
2330       for (auto& i_sigmaX : sector.m_sigmaX) {
2331         for (auto& iHist : i_sigmaX.second) {
2332           if (i_sigmaX.first == "sigmaXHit")
2333             iHist->Fill(hit.errXHit * 10000.);
2334           else if (i_sigmaX.first == "sigmaXTrk")
2335             iHist->Fill(hit.errXTrk * 10000.);
2336           else if (i_sigmaX.first == "sigmaX")
2337             iHist->Fill(hit.errX * 10000.);
2338         }
2339       }
2340       for (auto& i_sigmaY : sector.m_sigmaY) {
2341         for (auto& iHist : i_sigmaY.second) {
2342           if (i_sigmaY.first == "sigmaYHit")
2343             iHist->Fill(hit.errYHit * 10000.);
2344           else if (i_sigmaY.first == "sigmaYTrk")
2345             iHist->Fill(hit.errYTrk * 10000.);
2346           else if (i_sigmaY.first == "sigmaY")
2347             iHist->Fill(hit.errY * 10000.);
2348         }
2349       }
2350     }
2351   }
2352 }
2353 
2354 void ApeEstimator::fillHitHistsXForAnalyzerMode(const TrackStruct::HitParameterStruct& hit,
2355                                                 TrackerSectorStruct& sector) {
2356   std::map<std::string, TrackerSectorStruct::CorrelationHists>& m_corrHists(sector.m_correlationHistsX);
2357 
2358   // Cluster Parameters
2359   m_corrHists["WidthX"].fillCorrHistsX(hit, hit.widthX);
2360   m_corrHists["BaryStripX"].fillCorrHistsX(hit, hit.baryStripX);
2361 
2362   if (hit.isPixelHit) {
2363     m_corrHists["ChargePixel"].fillCorrHistsX(hit, hit.chargePixel);
2364     m_corrHists["ClusterProbXY"].fillCorrHistsX(hit, hit.clusterProbabilityXY);
2365     m_corrHists["ClusterProbQ"].fillCorrHistsX(hit, hit.clusterProbabilityQ);
2366     m_corrHists["ClusterProbXYQ"].fillCorrHistsX(hit, hit.clusterProbabilityXYQ);
2367     m_corrHists["LogClusterProb"].fillCorrHistsX(hit, hit.logClusterProbability);
2368     m_corrHists["IsOnEdge"].fillCorrHistsX(hit, hit.isOnEdge);
2369     m_corrHists["HasBadPixels"].fillCorrHistsX(hit, hit.hasBadPixels);
2370     m_corrHists["SpansTwoRoc"].fillCorrHistsX(hit, hit.spansTwoRoc);
2371     m_corrHists["QBin"].fillCorrHistsX(hit, hit.qBin);
2372 
2373   } else {
2374     m_corrHists["ChargeStrip"].fillCorrHistsX(hit, hit.chargeStrip);
2375     m_corrHists["MaxStrip"].fillCorrHistsX(hit, hit.maxStrip);
2376     m_corrHists["MaxCharge"].fillCorrHistsX(hit, hit.maxCharge);
2377     m_corrHists["MaxIndex"].fillCorrHistsX(hit, hit.maxIndex);
2378     m_corrHists["ChargeOnEdges"].fillCorrHistsX(hit, hit.chargeOnEdges);
2379     m_corrHists["ChargeAsymmetry"].fillCorrHistsX(hit, hit.chargeAsymmetry);
2380     m_corrHists["ChargeLRplus"].fillCorrHistsX(hit, hit.chargeLRplus);
2381     m_corrHists["ChargeLRminus"].fillCorrHistsX(hit, hit.chargeLRminus);
2382     m_corrHists["SOverN"].fillCorrHistsX(hit, hit.sOverN);
2383     m_corrHists["WidthProj"].fillCorrHistsX(hit, hit.projWidth);
2384     m_corrHists["WidthDiff"].fillCorrHistsX(hit, hit.projWidth - static_cast<float>(hit.widthX));
2385 
2386     sector.WidthVsWidthProjected->Fill(hit.projWidth, hit.widthX);
2387     sector.PWidthVsWidthProjected->Fill(hit.projWidth, hit.widthX);
2388 
2389     sector.WidthDiffVsMaxStrip->Fill(hit.maxStrip, hit.projWidth - static_cast<float>(hit.widthX));
2390     sector.PWidthDiffVsMaxStrip->Fill(hit.maxStrip, hit.projWidth - static_cast<float>(hit.widthX));
2391 
2392     sector.WidthDiffVsSigmaXHit->Fill(hit.errXHit, hit.projWidth - static_cast<float>(hit.widthX));
2393     sector.PWidthDiffVsSigmaXHit->Fill(hit.errXHit, hit.projWidth - static_cast<float>(hit.widthX));
2394 
2395     sector.WidthVsPhiSensX->Fill(hit.phiSensX * 180. / M_PI, hit.widthX);
2396     sector.PWidthVsPhiSensX->Fill(hit.phiSensX * 180. / M_PI, hit.widthX);
2397   }
2398 
2399   // Hit Parameters
2400   m_corrHists["SigmaXHit"].fillCorrHistsX(hit, hit.errXHit * 10000.);
2401   m_corrHists["SigmaXTrk"].fillCorrHistsX(hit, hit.errXTrk * 10000.);
2402   m_corrHists["SigmaX"].fillCorrHistsX(hit, hit.errX * 10000.);
2403 
2404   m_corrHists["PhiSens"].fillCorrHistsX(hit, hit.phiSens * 180. / M_PI);
2405   m_corrHists["PhiSensX"].fillCorrHistsX(hit, hit.phiSensX * 180. / M_PI);
2406   m_corrHists["PhiSensY"].fillCorrHistsX(hit, hit.phiSensY * 180. / M_PI);
2407 
2408   sector.XHit->Fill(hit.xHit);
2409   sector.XTrk->Fill(hit.xTrk);
2410   sector.SigmaX2->Fill(hit.errX2 * 10000. * 10000.);
2411 
2412   sector.ResX->Fill(hit.resX * 10000.);
2413   sector.NorResX->Fill(hit.norResX);
2414 
2415   sector.ProbX->Fill(hit.probX);
2416 
2417   sector.PhiSensXVsBarycentreX->Fill(hit.baryStripX, hit.phiSensX * 180. / M_PI);
2418   sector.PPhiSensXVsBarycentreX->Fill(hit.baryStripX, hit.phiSensX * 180. / M_PI);
2419 }
2420 
2421 void ApeEstimator::fillHitHistsYForAnalyzerMode(const TrackStruct::HitParameterStruct& hit,
2422                                                 TrackerSectorStruct& sector) {
2423   std::map<std::string, TrackerSectorStruct::CorrelationHists>& m_corrHists(sector.m_correlationHistsY);
2424   // Do not fill anything for strip
2425   if (!hit.isPixelHit)
2426     return;
2427 
2428   // Cluster Parameters
2429   m_corrHists["WidthY"].fillCorrHistsY(hit, hit.widthY);
2430   m_corrHists["BaryStripY"].fillCorrHistsY(hit, hit.baryStripY);
2431 
2432   m_corrHists["ChargePixel"].fillCorrHistsY(hit, hit.chargePixel);
2433   m_corrHists["ClusterProbXY"].fillCorrHistsY(hit, hit.clusterProbabilityXY);
2434   m_corrHists["ClusterProbQ"].fillCorrHistsY(hit, hit.clusterProbabilityQ);
2435   m_corrHists["ClusterProbXYQ"].fillCorrHistsY(hit, hit.clusterProbabilityXYQ);
2436   m_corrHists["LogClusterProb"].fillCorrHistsY(hit, hit.logClusterProbability);
2437   m_corrHists["IsOnEdge"].fillCorrHistsY(hit, hit.isOnEdge);
2438   m_corrHists["HasBadPixels"].fillCorrHistsY(hit, hit.hasBadPixels);
2439   m_corrHists["SpansTwoRoc"].fillCorrHistsY(hit, hit.spansTwoRoc);
2440   m_corrHists["QBin"].fillCorrHistsY(hit, hit.qBin);
2441 
2442   // Hit Parameters
2443   m_corrHists["SigmaYHit"].fillCorrHistsY(hit, hit.errYHit * 10000.);
2444   m_corrHists["SigmaYTrk"].fillCorrHistsY(hit, hit.errYTrk * 10000.);
2445   m_corrHists["SigmaY"].fillCorrHistsY(hit, hit.errY * 10000.);
2446 
2447   m_corrHists["PhiSens"].fillCorrHistsY(hit, hit.phiSens * 180. / M_PI);
2448   m_corrHists["PhiSensX"].fillCorrHistsY(hit, hit.phiSensX * 180. / M_PI);
2449   m_corrHists["PhiSensY"].fillCorrHistsY(hit, hit.phiSensY * 180. / M_PI);
2450 
2451   sector.YHit->Fill(hit.yHit);
2452   sector.YTrk->Fill(hit.yTrk);
2453   sector.SigmaY2->Fill(hit.errY2 * 10000. * 10000.);
2454 
2455   sector.ResY->Fill(hit.resY * 10000.);
2456   sector.NorResY->Fill(hit.norResY);
2457 
2458   sector.ProbY->Fill(hit.probY);
2459 
2460   sector.PhiSensYVsBarycentreY->Fill(hit.baryStripY, hit.phiSensY * 180. / M_PI);
2461   sector.PPhiSensYVsBarycentreY->Fill(hit.baryStripY, hit.phiSensY * 180. / M_PI);
2462 }
2463 
2464 void ApeEstimator::fillHistsForApeCalculation(const TrackStruct& trackStruct) {
2465   unsigned int goodHitsPerTrack(trackStruct.v_hitParams.size());
2466 
2467   if (parameterSet_.getParameter<bool>("applyTrackCuts")) {
2468     // which tracks to take? need min. nr. of selected hits?
2469     if (goodHitsPerTrack < minGoodHitsPerTrack_)
2470       return;
2471   }
2472 
2473   for (auto const& i_hit : trackStruct.v_hitParams) {
2474     // Put here from earlier method
2475     if (i_hit.hitState == TrackStruct::notAssignedToSectors)
2476       continue;
2477 
2478     for (auto& i_sector : m_tkSector_) {
2479       bool moduleInSector(false);
2480       for (auto const& i_hitSector : i_hit.v_sector) {
2481         if (i_sector.first == i_hitSector) {
2482           moduleInSector = true;
2483           break;
2484         }
2485       }
2486       if (!moduleInSector)
2487         continue;
2488 
2489       if (!calculateApe_)
2490         continue;
2491 
2492       if (i_hit.goodXMeasurement) {
2493         for (auto const& i_errBins : m_resErrBins_) {
2494           // Separate the bins for residual resolution w/o APE, to be consistent within iterations where APE will change (have same hit always in same bin)
2495           // So also fill this value in the histogram sigmaX
2496           // But of course use the normalized residual regarding the APE to have its influence in its width
2497           if (i_hit.errXWoApe < i_errBins.second.first || i_hit.errXWoApe >= i_errBins.second.second) {
2498             continue;
2499           }
2500           i_sector.second.m_binnedHists[i_errBins.first]["sigmaX"]->Fill(i_hit.errXWoApe);
2501           i_sector.second.m_binnedHists[i_errBins.first]["norResX"]->Fill(i_hit.norResX);
2502           break;
2503         }
2504         i_sector.second.ResX->Fill(i_hit.resX * 10000.);
2505         i_sector.second.NorResX->Fill(i_hit.norResX);
2506       }
2507 
2508       if (i_hit.goodYMeasurement) {
2509         for (auto const& i_errBins : m_resErrBins_) {
2510           // Separate the bins for residual resolution w/o APE, to be consistent within iterations where APE will change (have same hit always in same bin)
2511           // So also fill this value in the histogram sigmaY
2512           // But of course use the normalized residual regarding the APE to have its influence in its width
2513           if (i_hit.errYWoApe < i_errBins.second.first || i_hit.errYWoApe >= i_errBins.second.second) {
2514             continue;
2515           }
2516           i_sector.second.m_binnedHists[i_errBins.first]["sigmaY"]->Fill(i_hit.errYWoApe);
2517           i_sector.second.m_binnedHists[i_errBins.first]["norResY"]->Fill(i_hit.norResY);
2518           break;
2519         }
2520         i_sector.second.ResY->Fill(i_hit.resY * 10000.);
2521         i_sector.second.NorResY->Fill(i_hit.norResY);
2522       }
2523     }
2524   }
2525 }
2526 
2527 // -----------------------------------------------------------------------------------------------------------
2528 
2529 void ApeEstimator::calculateAPE() {
2530   // Loop over sectors for calculating APE
2531   for (auto& i_sector : m_tkSector_) {
2532     // Loop over residual error bins to calculate APE for every bin
2533     for (auto const& i_errBins : i_sector.second.m_binnedHists) {
2534       std::map<std::string, TH1*> m_Hists = i_errBins.second;
2535 
2536       // Fitting Parameters
2537       double integralX = m_Hists["norResX"]->Integral();
2538       i_sector.second.EntriesX->SetBinContent(i_errBins.first, integralX);
2539 
2540       if (i_sector.second.isPixel) {
2541         double integralY = m_Hists["norResY"]->Integral();
2542         i_sector.second.EntriesY->SetBinContent(i_errBins.first, integralY);
2543       }
2544     }
2545   }
2546 }
2547 
2548 // -----------------------------------------------------------------------------------------------------------
2549 
2550 bool ApeEstimator::isHit2D(const TrackingRecHit& hit) const {
2551   // we count SiStrip stereo modules as 2D if selected via countStereoHitAs2D_
2552   // (since they provide theta information)
2553   // --- NO, here it is always set to true ---
2554   if (!hit.isValid() || (hit.dimension() < 2 && !dynamic_cast<const SiStripRecHit1D*>(&hit))) {
2555     return false;  // real RecHit1D - but SiStripRecHit1D depends on countStereoHitAs2D_
2556   } else {
2557     const DetId detId(hit.geographicalId());
2558     if (detId.det() == DetId::Tracker) {
2559       if (detId.subdetId() == PixelSubdetector::PixelBarrel || detId.subdetId() == PixelSubdetector::PixelEndcap) {
2560         return true;  // pixel is always 2D
2561       } else {        // should be SiStrip now
2562         const SiStripDetId stripId(detId);
2563         if (stripId.stereo())
2564           return true;  // stereo modules
2565         else if (dynamic_cast<const SiStripRecHit1D*>(&hit) || dynamic_cast<const SiStripRecHit2D*>(&hit))
2566           return false;  // rphi modules hit
2567         //the following two are not used any more since ages...
2568         else if (dynamic_cast<const SiStripMatchedRecHit2D*>(&hit))
2569           return true;  // matched is 2D
2570         else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(&hit)) {
2571           const ProjectedSiStripRecHit2D* pH = static_cast<const ProjectedSiStripRecHit2D*>(&hit);
2572           return (this->isHit2D(pH->originalHit()));  // depends on original...
2573         } else {
2574           edm::LogError("UnknownType") << "@SUB=ApeEstimator::isHit2D"
2575                                        << "Tracker hit not in pixel, neither SiStripRecHit[12]D nor "
2576                                        << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
2577           return false;
2578         }
2579       }
2580     } else {  // not tracker??
2581       edm::LogWarning("DetectorMismatch") << "@SUB=ApeEstimator::isHit2D"
2582                                           << "Hit not in tracker with 'official' dimension >=2.";
2583       return true;  // dimension() >= 2 so accept that...
2584     }
2585   }
2586   // never reached...
2587 }
2588 
2589 // -----------------------------------------------------------------------------------------------------------
2590 
2591 // ------------ method called to for each event  ------------
2592 void ApeEstimator::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
2593   reco::BeamSpot beamSpot;
2594   edm::Handle<reco::BeamSpot> beamSpotHandle;
2595   iEvent.getByToken(offlinebeamSpot_, beamSpotHandle);
2596 
2597   if (beamSpotHandle.isValid()) {
2598     beamSpot = *beamSpotHandle;
2599   } else {
2600     edm::LogError("ApeEstimator") << "No beam spot available from EventSetup"
2601                                   << "\n...skip event";
2602     return;
2603   }
2604 
2605   edm::Handle<TrajTrackAssociationCollection> m_TrajTracksMap;
2606   iEvent.getByToken(tjTagToken_, m_TrajTracksMap);
2607 
2608   siStripClusterInfo_.initEvent(iSetup);
2609 
2610   if (analyzerMode_)
2611     tkDetector_.TrkSize->Fill(m_TrajTracksMap->size());
2612 
2613   if (maxTracksPerEvent_ != 0 && m_TrajTracksMap->size() > maxTracksPerEvent_)
2614     return;
2615 
2616   //Creation of (traj,track)
2617   typedef std::pair<const Trajectory*, const reco::Track*> ConstTrajTrackPair;
2618   typedef std::vector<ConstTrajTrackPair> ConstTrajTrackPairCollection;
2619   ConstTrajTrackPairCollection trajTracks;
2620 
2621   TrajTrackAssociationCollection::const_iterator i_trajTrack;
2622   for (i_trajTrack = m_TrajTracksMap->begin(); i_trajTrack != m_TrajTracksMap->end(); ++i_trajTrack) {
2623     trajTracks.push_back(ConstTrajTrackPair(&(*(*i_trajTrack).key), &(*(*i_trajTrack).val)));
2624   }
2625 
2626   //Loop over Tracks & Hits
2627   unsigned int trackSizeGood(0);
2628   ConstTrajTrackPairCollection::const_iterator iTrack;
2629   for (iTrack = trajTracks.begin(); iTrack != trajTracks.end(); ++iTrack) {
2630     const Trajectory* traj = (*iTrack).first;
2631     const reco::Track* track = (*iTrack).second;
2632 
2633     TrackStruct trackStruct;
2634     trackStruct.trkParams = this->fillTrackVariables(*track, *traj, beamSpot);
2635 
2636     if (trackCut_)
2637       continue;
2638 
2639     const std::vector<TrajectoryMeasurement> v_meas = (*traj).measurements();
2640 
2641     //Loop over Hits
2642     for (std::vector<TrajectoryMeasurement>::const_iterator i_meas = v_meas.begin(); i_meas != v_meas.end(); ++i_meas) {
2643       TrackStruct::HitParameterStruct hitParams = this->fillHitVariables(*i_meas, iSetup);
2644       if (this->hitSelected(hitParams))
2645         trackStruct.v_hitParams.push_back(hitParams);
2646     }
2647 
2648     if (analyzerMode_)
2649       this->fillHistsForAnalyzerMode(trackStruct);
2650     if (calculateApe_)
2651       this->fillHistsForApeCalculation(trackStruct);
2652 
2653     if (!trackStruct.v_hitParams.empty())
2654       ++trackSizeGood;
2655   }
2656   if (analyzerMode_ && trackSizeGood > 0)
2657     tkDetector_.TrkSizeGood->Fill(trackSizeGood);
2658 }
2659 
2660 // ------------ method called once each job just before starting event loop  ------------
2661 void ApeEstimator::beginJob() {
2662   this->hitSelection();
2663 
2664   this->sectorBuilder();
2665 
2666   this->residualErrorBinning();
2667 
2668   if (analyzerMode_)
2669     this->bookSectorHistsForAnalyzerMode();
2670 
2671   if (calculateApe_)
2672     this->bookSectorHistsForApeCalculation();
2673 
2674   if (analyzerMode_)
2675     this->bookTrackHists();
2676 }
2677 
2678 // ------------ method called once each job just after ending the event loop  ------------
2679 void ApeEstimator::endJob() {
2680   if (calculateApe_)
2681     this->calculateAPE();
2682 
2683   edm::LogInfo("HitSelector") << "\nThere are " << counter1 << " negative Errors calculated\n";
2684 }
2685 
2686 //define this as a plug-in
2687 DEFINE_FWK_MODULE(ApeEstimator);