File indexing completed on 2024-09-07 04:34:25
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022 #include <memory>
0023 #include <sstream>
0024 #include <fstream>
0025
0026
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
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
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
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
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
0213
0214
0215
0216
0217
0218
0219
0220
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
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
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
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 ;
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";
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.;
0631 double dzMax = zoomHists ? 15. : 100.;
0632 double pMax = zoomHists ? 200. : 2000.;
0633 double invPMax = zoomHists ? 0.05 : 10.;
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
0645 i_sector.second.Name = secDir.make<TH1F>("z_name", i_sector.second.name.c_str(), 1, 0, 1);
0646
0647
0648 if (i_sector.second.v_rawId.empty()) {
0649 continue;
0650 }
0651
0652 i_sector.second.setCorrHistParams(&secDir, norResXAbsMax, sigmaXHitMax, sigmaXMax);
0653
0654
0655 const bool pixelSector(i_sector.second.isPixel);
0656
0657
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
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.);
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.);
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
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
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
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
1098 i_sector.second.Name = secDir.make<TH1F>("z_name", i_sector.second.name.c_str(), 1, 0, 1);
1099
1100
1101 if (i_sector.second.v_rawId.empty()) {
1102 continue;
1103 }
1104
1105
1106 if (m_resErrBins_
1107 .empty()) {
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
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
1139
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
1152
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.;
1177 double dzmax = zoomHists ? 15. : 100.;
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();
1285 trkParams.hitsInvalid = trkParams.hitsSize - trkParams.hitsValid;
1286 trkParams.layersMissed =
1287 track.lost();
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 ||
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);
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
1391 const LocalPoint& lPHit = hit.localPosition();
1392 const LocalPoint& lPTrk = tsos.localPosition();
1393
1394
1395
1396
1397 const LocalError& errHitApe = hit.localPositionError();
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 }
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
1466
1467
1468
1469
1470
1471
1472 hitParams.hitState = TrackStruct::negativeError;
1473 return hitParams;
1474 }
1475
1476
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
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
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
1655 if (!hit.detUnit()) {
1656 hitParams.hitState = TrackStruct::invalid;
1657 return hitParams;
1658 }
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.;
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
1738
1739
1740
1741
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
1749 vPE2 = std::make_pair(TrackStruct::ok, this->rectangularPositionAndError2(localPoint, localError));
1750 } else if (subdetId == StripSubdetector::TID || subdetId == StripSubdetector::TEC) {
1751
1752 if (!hit.detUnit())
1753 return vPE2;
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
1763
1764
1765
1766
1767
1768
1769
1770
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
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
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();
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();
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
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
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 {
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
2114
2115
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 {
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
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
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
2283 this->fillHitHistsXForAnalyzerMode(hit, sector);
2284
2285
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
2308 this->fillHitHistsYForAnalyzerMode(hit, sector);
2309
2310
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
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
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
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
2425 if (!hit.isPixelHit)
2426 return;
2427
2428
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
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
2469 if (goodHitsPerTrack < minGoodHitsPerTrack_)
2470 return;
2471 }
2472
2473 for (auto const& i_hit : trackStruct.v_hitParams) {
2474
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
2495
2496
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
2511
2512
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
2531 for (auto& i_sector : m_tkSector_) {
2532
2533 for (auto const& i_errBins : i_sector.second.m_binnedHists) {
2534 std::map<std::string, TH1*> m_Hists = i_errBins.second;
2535
2536
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
2552
2553
2554 if (!hit.isValid() || (hit.dimension() < 2 && !dynamic_cast<const SiStripRecHit1D*>(&hit))) {
2555 return false;
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;
2561 } else {
2562 const SiStripDetId stripId(detId);
2563 if (stripId.stereo())
2564 return true;
2565 else if (dynamic_cast<const SiStripRecHit1D*>(&hit) || dynamic_cast<const SiStripRecHit2D*>(&hit))
2566 return false;
2567
2568 else if (dynamic_cast<const SiStripMatchedRecHit2D*>(&hit))
2569 return true;
2570 else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(&hit)) {
2571 const ProjectedSiStripRecHit2D* pH = static_cast<const ProjectedSiStripRecHit2D*>(&hit);
2572 return (this->isHit2D(pH->originalHit()));
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 {
2581 edm::LogWarning("DetectorMismatch") << "@SUB=ApeEstimator::isHit2D"
2582 << "Hit not in tracker with 'official' dimension >=2.";
2583 return true;
2584 }
2585 }
2586
2587 }
2588
2589
2590
2591
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
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
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
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
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
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
2687 DEFINE_FWK_MODULE(ApeEstimator);
