File indexing completed on 2024-04-06 12:26:18
0001
0002 #include <map>
0003 #include <vector>
0004 #include <algorithm>
0005
0006
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/Utilities/interface/InputTag.h"
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0018 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0019 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0020 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0021 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0022 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0023 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0024 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0025 #include "DataFormats/Common/interface/DetSetVector.h"
0026 #include "DataFormats/DetId/interface/DetId.h"
0027 #include "DataFormats/Common/interface/Handle.h"
0028 #include "DataFormats/Phase2TrackerCluster/interface/Phase2TrackerCluster1D.h"
0029 #include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h"
0030 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
0031 #include "RecoLocalTracker/Phase2TrackerRecHits/interface/Phase2StripCPE.h"
0032 #include "SimDataFormats/TrackerDigiSimLink/interface/PixelDigiSimLink.h"
0033 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0034 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0035 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0036 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0037 #include "CommonTools/Utils/interface/TFileDirectory.h"
0038 #include "DataFormats/TrackerRecHit2D/interface/VectorHit.h"
0039 #include "DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h"
0040 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0041 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0042 #include "SimTracker/TrackerHitAssociation/interface/ClusterTPAssociation.h"
0043 #include "Geometry/CommonDetUnit/interface/StackGeomDet.h"
0044 #include "DataFormats/Phase2TrackerDigi/interface/Phase2TrackerDigi.h"
0045 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0046 #include "RecoLocalTracker/Records/interface/TkPhase2OTCPERecord.h"
0047 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0048 #include "RecoLocalTracker/ClusterParameterEstimator/interface/ClusterParameterEstimator.h"
0049
0050
0051 #include <TH1F.h>
0052 #include <TH2D.h>
0053 #include <TGraph.h>
0054 #include <THStack.h>
0055 #include <TCanvas.h>
0056 #include <TTree.h>
0057 #include <TArrow.h>
0058
0059 struct VHHistos {
0060 THStack* numberVHsMixed;
0061 TH1F* numberVHsPS;
0062 TH1F* numberVHs2S;
0063
0064 TGraph* globalPosXY[3];
0065 TGraph* localPosXY[3];
0066
0067 TH1F* deltaXVHSimHits[3];
0068 TH1F* deltaYVHSimHits[3];
0069
0070 TH1F* deltaXVHSimHits_P[3];
0071 TH1F* deltaYVHSimHits_P[3];
0072
0073 TH1F* digiEfficiency[3];
0074
0075 TH1F* totalSimHits;
0076 TH1F* primarySimHits;
0077 TH1F* otherSimHits;
0078
0079 TH1F* curvature;
0080 TH1F* width;
0081 TH1F* deltaXlocal;
0082 };
0083
0084 class VectorHitsBuilderValidation : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0085 public:
0086 typedef edm::Ref<edmNew::DetSetVector<Phase2TrackerCluster1D>, Phase2TrackerCluster1D> Phase2TrackerCluster1DRef;
0087
0088 typedef std::map<unsigned int, std::vector<PSimHit> > SimHitsMap;
0089 typedef std::map<unsigned int, SimTrack> SimTracksMap;
0090
0091 explicit VectorHitsBuilderValidation(const edm::ParameterSet&);
0092 ~VectorHitsBuilderValidation();
0093 void beginJob();
0094 void analyze(const edm::Event&, const edm::EventSetup&);
0095
0096 static void fillDescriptions(edm::ConfigurationDescriptions&);
0097
0098 private:
0099 std::map<unsigned int, VHHistos>::iterator createLayerHistograms(unsigned int);
0100 void CreateVHsXYGraph(const std::vector<Global3DPoint>, const std::vector<Global3DVector>);
0101 void CreateVHsRZGraph(const std::vector<Global3DPoint>, const std::vector<Global3DVector>);
0102 void CreateWindowCorrGraph();
0103
0104 unsigned int getLayerNumber(const DetId&);
0105 unsigned int getModuleNumber(const DetId& detid);
0106 void printCluster(const GeomDetUnit* geomDetUnit, const OmniClusterRef cluster);
0107
0108 std::pair<bool, uint32_t> isTrue(const VectorHit vh,
0109 const edm::Handle<edm::DetSetVector<PixelDigiSimLink> >& siphase2SimLinks,
0110 DetId& detId) const;
0111 std::vector<std::pair<uint32_t, EncodedEventId> > getSimTrackIds(
0112 const edm::Handle<edm::DetSetVector<PixelDigiSimLink> >&, const DetId&, uint32_t) const;
0113 unsigned int getSimTrackId(const edm::Handle<edm::DetSetVector<PixelDigiSimLink> >& pixelSimLinks,
0114 const DetId& detId,
0115 unsigned int channel) const;
0116
0117 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomEsToken_;
0118 const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoEsToken_;
0119 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldEsToken_;
0120 const edm::ESGetToken<ClusterParameterEstimator<Phase2TrackerCluster1D>, TkPhase2OTCPERecord> cpeEsToken_;
0121
0122 edm::EDGetTokenT<edmNew::DetSetVector<Phase2TrackerCluster1D> > srcClu_;
0123 edm::EDGetTokenT<VectorHitCollection> VHacc_;
0124 edm::EDGetTokenT<VectorHitCollection> VHrej_;
0125 edm::EDGetTokenT<edm::DetSetVector<PixelDigiSimLink> > siphase2OTSimLinksToken_;
0126 edm::EDGetTokenT<edm::PSimHitContainer> simHitsToken_;
0127 edm::EDGetTokenT<edm::SimTrackContainer> simTracksToken_;
0128 edm::EDGetTokenT<edm::SimVertexContainer> simVerticesToken_;
0129 edm::EDGetTokenT<TrackingParticleCollection> trackingParticleToken_;
0130
0131 const TrackerGeometry* tkGeom_;
0132 const TrackerTopology* tkTopo_;
0133 const MagneticField* magField_;
0134 const ClusterParameterEstimator<Phase2TrackerCluster1D>* cpe_;
0135
0136 TTree* tree_;
0137 TGraph* trackerLayoutRZ_[3];
0138 TGraph* trackerLayoutXY_[3];
0139 TGraph* trackerLayoutXYBar_;
0140 TGraph* trackerLayoutXYEC_;
0141 TGraph* localPosXvsDeltaX_[3];
0142 TGraph* localPosYvsDeltaY_[3];
0143 TCanvas* VHXY_;
0144 TCanvas* VHRZ_;
0145 std::vector<TArrow*> arrowVHs_;
0146
0147 TH2D* ParallaxCorrectionRZ_;
0148 TH1F* VHaccLayer_;
0149 TH1F* VHrejLayer_;
0150 TH1F* VHaccTrueLayer_;
0151 TH1F* VHrejTrueLayer_;
0152 TH1F* VHaccTrue_signal_Layer_;
0153 TH1F* VHrejTrue_signal_Layer_;
0154
0155 std::map<unsigned int, VHHistos> histograms_;
0156 };
0157
0158 VectorHitsBuilderValidation::VectorHitsBuilderValidation(const edm::ParameterSet& conf)
0159 : geomEsToken_(esConsumes()),
0160 topoEsToken_(esConsumes()),
0161 magFieldEsToken_(esConsumes()),
0162 cpeEsToken_(esConsumes(conf.getParameter<edm::ESInputTag>("CPE"))) {
0163 srcClu_ =
0164 consumes<edmNew::DetSetVector<Phase2TrackerCluster1D> >(edm::InputTag(conf.getParameter<std::string>("src")));
0165 VHacc_ = consumes<VectorHitCollection>(edm::InputTag(conf.getParameter<edm::InputTag>("VH_acc")));
0166 VHrej_ = consumes<VectorHitCollection>(edm::InputTag(conf.getParameter<edm::InputTag>("VH_rej")));
0167 siphase2OTSimLinksToken_ = consumes<edm::DetSetVector<PixelDigiSimLink> >(conf.getParameter<edm::InputTag>("links"));
0168 simHitsToken_ = consumes<edm::PSimHitContainer>(edm::InputTag("g4SimHits", "TrackerHitsPixelBarrelLowTof"));
0169 simTracksToken_ = consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"));
0170 simVerticesToken_ = consumes<edm::SimVertexContainer>(edm::InputTag("g4SimHits"));
0171 trackingParticleToken_ =
0172 consumes<TrackingParticleCollection>(conf.getParameter<edm::InputTag>("trackingParticleSrc"));
0173 }
0174
0175 VectorHitsBuilderValidation::~VectorHitsBuilderValidation() = default;
0176
0177 void VectorHitsBuilderValidation::beginJob() {
0178 edm::Service<TFileService> fs;
0179 fs->file().cd("/");
0180 TFileDirectory td = fs->mkdir("Common");
0181
0182
0183 tree_ = td.make<TTree>("VectorHits", "VectorHits");
0184
0185
0186 TFileDirectory tdGloPos = td.mkdir("GlobalPositions");
0187 trackerLayoutRZ_[0] = tdGloPos.make<TGraph>();
0188 trackerLayoutRZ_[0]->SetName("RVsZ_Mixed");
0189 trackerLayoutRZ_[1] = tdGloPos.make<TGraph>();
0190 trackerLayoutRZ_[1]->SetName("RVsZ_Pixel");
0191 trackerLayoutRZ_[2] = tdGloPos.make<TGraph>();
0192 trackerLayoutRZ_[2]->SetName("RVsZ_Strip");
0193 trackerLayoutXY_[0] = tdGloPos.make<TGraph>();
0194 trackerLayoutXY_[0]->SetName("YVsX_Mixed");
0195 trackerLayoutXY_[1] = tdGloPos.make<TGraph>();
0196 trackerLayoutXY_[1]->SetName("YVsX_Pixel");
0197 trackerLayoutXY_[2] = tdGloPos.make<TGraph>();
0198 trackerLayoutXY_[2]->SetName("YVsX_Strip");
0199 trackerLayoutXYBar_ = tdGloPos.make<TGraph>();
0200 trackerLayoutXYBar_->SetName("YVsXBar");
0201 trackerLayoutXYEC_ = tdGloPos.make<TGraph>();
0202 trackerLayoutXYEC_->SetName("YVsXEC");
0203
0204 TFileDirectory tdLocPos = td.mkdir("LocalPositions");
0205 localPosXvsDeltaX_[0] = tdLocPos.make<TGraph>();
0206 localPosXvsDeltaX_[0]->SetName("localPosXvsDeltaX_Mixed");
0207 localPosXvsDeltaX_[1] = tdLocPos.make<TGraph>();
0208 localPosXvsDeltaX_[1]->SetName("localPosXvsDeltaX_Pixel");
0209 localPosXvsDeltaX_[2] = tdLocPos.make<TGraph>();
0210 localPosXvsDeltaX_[2]->SetName("localPosXvsDeltaX_Strip");
0211 localPosYvsDeltaY_[0] = tdLocPos.make<TGraph>();
0212 localPosYvsDeltaY_[0]->SetName("localPosYvsDeltaY_Mixed");
0213 localPosYvsDeltaY_[1] = tdLocPos.make<TGraph>();
0214 localPosYvsDeltaY_[1]->SetName("localPosYvsDeltaY_Pixel");
0215 localPosYvsDeltaY_[2] = tdLocPos.make<TGraph>();
0216 localPosYvsDeltaY_[2]->SetName("localPosYvsDeltaY_Strip");
0217
0218
0219 TFileDirectory tdArr = td.mkdir("Directions");
0220
0221 TFileDirectory tdWid = td.mkdir("CombinatorialStudies");
0222 ParallaxCorrectionRZ_ =
0223 tdWid.make<TH2D>("ParallaxCorrectionRZ", "ParallaxCorrectionRZ", 100, 0., 300., 100., 0., 120.);
0224 ParallaxCorrectionRZ_->SetName("ParallaxCorrectionFactor");
0225 VHaccLayer_ = tdWid.make<TH1F>("VHacceptedLayer", "VHacceptedLayer", 250, 0., 250.);
0226 VHaccLayer_->SetName("VHaccepted");
0227 VHrejLayer_ = tdWid.make<TH1F>("VHrejectedLayer", "VHrejectedLayer", 250, 0., 250.);
0228 VHrejLayer_->SetName("VHrejected");
0229 VHaccTrueLayer_ = tdWid.make<TH1F>("VHaccTrueLayer", "VHaccTrueLayer", 250, 0., 250.);
0230 VHaccTrueLayer_->SetName("VHaccepted_true");
0231 VHrejTrueLayer_ = tdWid.make<TH1F>("VHrejTrueLayer", "VHrejTrueLayer", 250, 0., 250.);
0232 VHrejTrueLayer_->SetName("VHrejected_true");
0233 VHaccTrue_signal_Layer_ = tdWid.make<TH1F>("VHaccTrueSignalLayer", "VHaccTrueSignalLayer", 250, 0., 250.);
0234 VHaccTrue_signal_Layer_->SetName("VHaccepted_true_signal");
0235 VHrejTrue_signal_Layer_ = tdWid.make<TH1F>("VHrejTrueSignalLayer", "VHrejTrueSignalLayer", 250, 0., 250.);
0236 VHrejTrue_signal_Layer_->SetName("VHrejected_true_signal");
0237 }
0238
0239 void VectorHitsBuilderValidation::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0240
0241
0242
0243 edm::Handle<Phase2TrackerCluster1DCollectionNew> clusters;
0244 event.getByToken(srcClu_, clusters);
0245
0246
0247 edm::Handle<VectorHitCollection> vhsAcc;
0248 event.getByToken(VHacc_, vhsAcc);
0249
0250 edm::Handle<VectorHitCollection> vhsRej;
0251 event.getByToken(VHrej_, vhsRej);
0252
0253
0254 cpe_ = &eventSetup.getData(cpeEsToken_);
0255
0256
0257 edm::Handle<edm::DetSetVector<PixelDigiSimLink> > siphase2SimLinks;
0258 event.getByToken(siphase2OTSimLinksToken_, siphase2SimLinks);
0259
0260
0261 edm::Handle<edm::PSimHitContainer> simHitsRaw;
0262 event.getByToken(simHitsToken_, simHitsRaw);
0263
0264
0265 edm::Handle<edm::SimTrackContainer> simTracksRaw;
0266 event.getByToken(simTracksToken_, simTracksRaw);
0267
0268
0269 edm::Handle<edm::SimVertexContainer> simVertices;
0270 event.getByToken(simVerticesToken_, simVertices);
0271
0272
0273 tkGeom_ = &eventSetup.getData(geomEsToken_);
0274
0275
0276 tkTopo_ = &eventSetup.getData(topoEsToken_);
0277
0278
0279 magField_ = &eventSetup.getData(magFieldEsToken_);
0280
0281
0282 edm::Handle<TrackingParticleCollection> TPCollectionH;
0283 event.getByToken(trackingParticleToken_, TPCollectionH);
0284
0285 auto clusterTPList = std::make_unique<ClusterTPAssociation>(TPCollectionH);
0286 std::map<std::pair<size_t, EncodedEventId>, TrackingParticleRef> mapping;
0287
0288 for (TrackingParticleCollection::size_type itp = 0; itp < TPCollectionH.product()->size(); ++itp) {
0289 TrackingParticleRef trackingParticle(TPCollectionH, itp);
0290 EncodedEventId eid(trackingParticle->eventId());
0291 for (std::vector<SimTrack>::const_iterator itrk = trackingParticle->g4Track_begin();
0292 itrk != trackingParticle->g4Track_end();
0293 ++itrk) {
0294 std::pair<uint32_t, EncodedEventId> trkid(itrk->trackId(), eid);
0295 LogTrace("VectorHitsBuilderValidation")
0296 << "creating map for id: " << trkid.first << " with tp: " << trackingParticle.key();
0297 mapping.insert(std::make_pair(trkid, trackingParticle));
0298 }
0299 }
0300
0301
0302 int eventNum;
0303 int layer;
0304 int module_id;
0305 int module_number;
0306 int module_type;
0307 int VHacc = 0.0;
0308 int VHrej = 0.0;
0309 int vh_isTrue;
0310
0311 float x_global, y_global, z_global;
0312 float vh_x_local, vh_y_local;
0313 float vh_x_le, vh_y_le;
0314 float curvature, phi;
0315 float QOverPT, QOverP;
0316 float chi2;
0317
0318 int low_tp_id, upp_tp_id;
0319 float vh_sim_trackPt;
0320 float sim_x_local, sim_y_local;
0321 float sim_x_global, sim_y_global, sim_z_global;
0322 float low_x_global, low_y_global, low_z_global;
0323 float upp_x_global, upp_y_global, upp_z_global;
0324 float low_xx_global_err, low_yy_global_err, low_zz_global_err;
0325 float low_xy_global_err, low_zx_global_err, low_zy_global_err;
0326 float upp_xx_global_err, upp_yy_global_err, upp_zz_global_err;
0327 float upp_xy_global_err, upp_zx_global_err, upp_zy_global_err;
0328 float deltaXVHSimHits, deltaYVHSimHits;
0329 int multiplicity;
0330 float width, deltaXlocal;
0331 unsigned int processType(99);
0332
0333 tree_->Branch("event", &eventNum, "eventNum/I");
0334 tree_->Branch("accepted", &VHacc, "VHacc/I");
0335 tree_->Branch("rejected", &VHrej, "VHrej/I");
0336 tree_->Branch("layer", &layer, "layer/I");
0337 tree_->Branch("module_id", &module_id, "module_id/I");
0338 tree_->Branch("module_type", &module_type, "module_type/I");
0339 tree_->Branch("module_number", &module_number, "module_number/I");
0340 tree_->Branch("vh_isTrue", &vh_isTrue, "vh_isTrue/I");
0341 tree_->Branch("x_global", &x_global, "x_global/F");
0342 tree_->Branch("y_global", &y_global, "y_global/F");
0343 tree_->Branch("z_global", &z_global, "z_global/F");
0344 tree_->Branch("vh_x_local", &vh_x_local, "vh_x_local/F");
0345 tree_->Branch("vh_y_local", &vh_y_local, "vh_y_local/F");
0346 tree_->Branch("vh_x_lError", &vh_x_le, "vh_x_le/F");
0347 tree_->Branch("vh_y_lError", &vh_y_le, "vh_y_le/F");
0348 tree_->Branch("curvature", &curvature, "curvature/F");
0349 tree_->Branch("chi2", &chi2, "chi2/F");
0350 tree_->Branch("phi", &phi, "phi/F");
0351 tree_->Branch("QOverP", &QOverP, "QOverP/F");
0352 tree_->Branch("QOverPT", &QOverPT, "QOverPT/F");
0353 tree_->Branch("low_tp_id", &low_tp_id, "low_tp_id/I");
0354 tree_->Branch("upp_tp_id", &upp_tp_id, "upp_tp_id/I");
0355 tree_->Branch("vh_sim_trackPt", &vh_sim_trackPt, "vh_sim_trackPt/F");
0356 tree_->Branch("sim_x_local", &sim_x_local, "sim_x_local/F");
0357 tree_->Branch("sim_y_local", &sim_y_local, "sim_y_local/F");
0358 tree_->Branch("sim_x_global", &sim_x_global, "sim_x_global/F");
0359 tree_->Branch("sim_y_global", &sim_y_global, "sim_y_global/F");
0360 tree_->Branch("sim_z_global", &sim_z_global, "sim_z_global/F");
0361 tree_->Branch("low_x_global", &low_x_global, "low_x_global/F");
0362 tree_->Branch("low_y_global", &low_y_global, "low_y_global/F");
0363 tree_->Branch("low_z_global", &low_z_global, "low_z_global/F");
0364 tree_->Branch("low_xx_global_err", &low_xx_global_err, "low_xx_global_err/F");
0365 tree_->Branch("low_yy_global_err", &low_yy_global_err, "low_yy_global_err/F");
0366 tree_->Branch("low_zz_global_err", &low_zz_global_err, "low_zz_global_err/F");
0367 tree_->Branch("low_xy_global_err", &low_xy_global_err, "low_xy_global_err/F");
0368 tree_->Branch("low_zx_global_err", &low_zx_global_err, "low_zx_global_err/F");
0369 tree_->Branch("low_zy_global_err", &low_zy_global_err, "low_zy_global_err/F");
0370 tree_->Branch("upp_x_global", &upp_x_global, "upp_x_global/F");
0371 tree_->Branch("upp_y_global", &upp_y_global, "upp_y_global/F");
0372 tree_->Branch("upp_z_global", &upp_z_global, "upp_z_global/F");
0373 tree_->Branch("upp_xx_global_err", &upp_xx_global_err, "upp_xx_global_err/F");
0374 tree_->Branch("upp_yy_global_err", &upp_yy_global_err, "upp_yy_global_err/F");
0375 tree_->Branch("upp_zz_global_err", &upp_zz_global_err, "upp_zz_global_err/F");
0376 tree_->Branch("upp_xy_global_err", &upp_xy_global_err, "upp_xy_global_err/F");
0377 tree_->Branch("upp_zx_global_err", &upp_zx_global_err, "upp_zx_global_err/F");
0378 tree_->Branch("upp_zy_global_err", &upp_zy_global_err, "upp_zy_global_err/F");
0379 tree_->Branch("deltaXVHSimHits", &deltaXVHSimHits, "deltaXVHSimHits/F");
0380 tree_->Branch("deltaYVHSimHits", &deltaYVHSimHits, "deltaYVHSimHits/F");
0381 tree_->Branch("multiplicity", &multiplicity, "multiplicity/I");
0382 tree_->Branch("width", &width, "width/F");
0383 tree_->Branch("deltaXlocal", &deltaXlocal, "deltaXlocal/F");
0384 tree_->Branch("processType", &processType, "processType/i");
0385
0386
0387 SimTracksMap simTracks;
0388 for (const auto& simTrackIt : *simTracksRaw.product())
0389 simTracks.emplace(std::pair<unsigned int, SimTrack>(simTrackIt.trackId(), simTrackIt));
0390
0391
0392 SimHitsMap simHitsDetUnit;
0393 SimHitsMap simHitsTrackId;
0394 for (const auto& simHitIt : *simHitsRaw.product()) {
0395 SimHitsMap::iterator simHitsDetUnitIt(simHitsDetUnit.find(simHitIt.detUnitId()));
0396 if (simHitsDetUnitIt == simHitsDetUnit.end()) {
0397 std::pair<SimHitsMap::iterator, bool> newIt(simHitsDetUnit.insert(
0398 std::pair<unsigned int, std::vector<PSimHit> >(simHitIt.detUnitId(), std::vector<PSimHit>())));
0399 simHitsDetUnitIt = newIt.first;
0400 }
0401 simHitsDetUnitIt->second.push_back(simHitIt);
0402
0403 SimHitsMap::iterator simHitsTrackIdIt(simHitsTrackId.find(simHitIt.trackId()));
0404 if (simHitsTrackIdIt == simHitsTrackId.end()) {
0405 std::pair<SimHitsMap::iterator, bool> newIt(simHitsTrackId.insert(
0406 std::pair<unsigned int, std::vector<PSimHit> >(simHitIt.trackId(), std::vector<PSimHit>())));
0407 simHitsTrackIdIt = newIt.first;
0408 }
0409 simHitsTrackIdIt->second.push_back(simHitIt);
0410 }
0411
0412
0413 for (const auto& DSViter : *clusters) {
0414 unsigned int rawid(DSViter.detId());
0415 DetId detId(rawid);
0416 const GeomDetUnit* geomDetUnit(tkGeom_->idToDetUnit(detId));
0417 const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>(geomDetUnit);
0418 for (const auto& clustIt : DSViter) {
0419 auto&& lparams = cpe_->localParameters(clustIt, *theGeomDet);
0420 Global3DPoint gparams = theGeomDet->surface().toGlobal(lparams.first);
0421 LogTrace("VectorHitsBuilderValidation") << "phase2 OT clusters: " << gparams << " DetId: " << rawid;
0422 }
0423 }
0424
0425 for (const auto& DSViter : *vhsAcc) {
0426 for (const auto& vhIt : DSViter) {
0427 LogTrace("VectorHitsBuilderValidation") << "accepted VH: " << vhIt;
0428 }
0429 }
0430 for (const auto& DSViter : *vhsRej) {
0431 for (const auto& vhIt : DSViter) {
0432 LogTrace("VectorHitsBuilderValidation") << "rejected VH: " << vhIt;
0433 }
0434 }
0435
0436 eventNum = event.id().event();
0437
0438 unsigned int nVHsTot(0), nVHsPSTot(0), nVHs2STot(0);
0439 std::vector<Global3DPoint> glVHs;
0440 std::vector<Global3DVector> dirVHs;
0441 std::vector<int> detIds;
0442
0443
0444 for (const auto DSViter : *vhsAcc) {
0445
0446 unsigned int rawid(DSViter.detId());
0447 module_id = rawid;
0448 DetId detId(rawid);
0449
0450 module_number = getModuleNumber(detId);
0451 layer = getLayerNumber(detId);
0452
0453 LogDebug("VectorHitsBuilderValidation") << "Layer: " << layer << " det id" << rawid << std::endl;
0454
0455
0456 const GeomDet* geomDet(tkGeom_->idToDet(detId));
0457 if (!geomDet)
0458 break;
0459
0460
0461 std::map<unsigned int, VHHistos>::iterator histogramLayer(histograms_.find(layer));
0462 if (histogramLayer == histograms_.end())
0463 histogramLayer = createLayerHistograms(layer);
0464
0465 unsigned int nVHsPS(0), nVHs2S(0);
0466
0467 LogDebug("VectorHitsBuilderValidation") << "DSViter size: " << DSViter.size();
0468
0469
0470 for (const auto& vhIt : DSViter) {
0471
0472 if (vhIt.isValid()) {
0473 LogDebug("VectorHitsBuilderValidation") << " vh analyzing ...";
0474 chi2 = vhIt.chi2();
0475 LogTrace("VectorHitsBuilderValidation") << "VH chi2 " << chi2 << std::endl;
0476
0477 Local3DPoint localPosVH = vhIt.localPosition();
0478 vh_x_local = localPosVH.x();
0479 vh_y_local = localPosVH.y();
0480 LogTrace("VectorHitsBuilderValidation") << "local VH position " << localPosVH << std::endl;
0481
0482 LocalError localErrVH = vhIt.localPositionError();
0483 vh_x_le = localErrVH.xx();
0484 vh_y_le = localErrVH.yy();
0485 LogTrace("VectorHitsBuilderValidation") << "local VH error " << localErrVH << std::endl;
0486
0487 Global3DPoint globalPosVH = geomDet->surface().toGlobal(localPosVH);
0488 x_global = globalPosVH.x();
0489 y_global = globalPosVH.y();
0490 z_global = globalPosVH.z();
0491 glVHs.push_back(globalPosVH);
0492 LogTrace("VectorHitsBuilderValidation") << " global VH position " << globalPosVH << std::endl;
0493
0494 Local3DVector localDirVH = vhIt.localDirection();
0495 LogTrace("VectorHitsBuilderValidation") << "local VH direction " << localDirVH << std::endl;
0496
0497 VectorHit vh = vhIt;
0498 Global3DVector globalDirVH = vh.globalDirectionVH();
0499 dirVHs.push_back(globalDirVH);
0500 LogTrace("VectorHitsBuilderValidation") << "global VH direction " << globalDirVH << std::endl;
0501
0502
0503 trackerLayoutRZ_[0]->SetPoint(nVHsTot, globalPosVH.z(), globalPosVH.perp());
0504 trackerLayoutXY_[0]->SetPoint(nVHsTot, globalPosVH.x(), globalPosVH.y());
0505
0506 if (layer < 100)
0507 trackerLayoutXYBar_->SetPoint(nVHsTot, globalPosVH.x(), globalPosVH.y());
0508 else
0509 trackerLayoutXYEC_->SetPoint(nVHsTot, globalPosVH.x(), globalPosVH.y());
0510
0511 histogramLayer->second.localPosXY[0]->SetPoint(nVHsTot, vh_x_local, vh_y_local);
0512 histogramLayer->second.globalPosXY[0]->SetPoint(nVHsTot, globalPosVH.x(), globalPosVH.y());
0513
0514 localPosXvsDeltaX_[0]->SetPoint(nVHsTot, vh_x_local, localDirVH.x());
0515 localPosYvsDeltaY_[0]->SetPoint(nVHsTot, vh_y_local, localDirVH.y());
0516
0517
0518 const StackGeomDet* stackDet = dynamic_cast<const StackGeomDet*>(geomDet);
0519 const PixelGeomDetUnit* geomDetLower = dynamic_cast<const PixelGeomDetUnit*>(stackDet->lowerDet());
0520 DetId lowerDetId = stackDet->lowerDet()->geographicalId();
0521 DetId upperDetId = stackDet->upperDet()->geographicalId();
0522
0523 TrackerGeometry::ModuleType mType = tkGeom_->getDetectorType(lowerDetId);
0524 module_type = 0;
0525 if (mType == TrackerGeometry::ModuleType::Ph2PSP) {
0526 module_type = 1;
0527 trackerLayoutRZ_[1]->SetPoint(nVHsPSTot, globalPosVH.z(), globalPosVH.perp());
0528 trackerLayoutXY_[1]->SetPoint(nVHsPSTot, globalPosVH.x(), globalPosVH.y());
0529
0530 histogramLayer->second.localPosXY[1]->SetPoint(nVHsPSTot, vh_x_local, vh_y_local);
0531 histogramLayer->second.globalPosXY[1]->SetPoint(nVHsPSTot, globalPosVH.x(), globalPosVH.y());
0532
0533 localPosXvsDeltaX_[1]->SetPoint(nVHsPSTot, vh_x_local, localDirVH.x());
0534 localPosYvsDeltaY_[1]->SetPoint(nVHsPSTot, vh_y_local, localDirVH.y());
0535
0536 ++nVHsPS;
0537 ++nVHsPSTot;
0538 }
0539
0540
0541 else if (mType == TrackerGeometry::ModuleType::Ph2SS) {
0542 module_type = 2;
0543 trackerLayoutRZ_[2]->SetPoint(nVHs2STot, globalPosVH.z(), globalPosVH.perp());
0544 trackerLayoutXY_[2]->SetPoint(nVHs2STot, globalPosVH.x(), globalPosVH.y());
0545
0546 histogramLayer->second.localPosXY[2]->SetPoint(nVHs2STot, vh_x_local, vh_y_local);
0547 histogramLayer->second.globalPosXY[2]->SetPoint(nVHs2STot, globalPosVH.x(), globalPosVH.y());
0548
0549 localPosXvsDeltaX_[2]->SetPoint(nVHs2STot, vh_x_local, localDirVH.x());
0550 localPosYvsDeltaY_[2]->SetPoint(nVHs2STot, vh_y_local, localDirVH.y());
0551
0552 ++nVHs2S;
0553 ++nVHs2STot;
0554 } else if (mType == TrackerGeometry::ModuleType::Ph2PSS) {
0555 edm::LogError("VectorHitsBuilderValidation") << "module type " << module_type << " should never happen!";
0556 }
0557 LogTrace("VectorHitsBuilderValidation") << "module type " << module_type << std::endl;
0558
0559
0560 low_x_global = vhIt.lowerGlobalPos().x();
0561 low_y_global = vhIt.lowerGlobalPos().y();
0562 low_z_global = vhIt.lowerGlobalPos().z();
0563 upp_x_global = vhIt.upperGlobalPos().x();
0564 upp_y_global = vhIt.upperGlobalPos().y();
0565 upp_z_global = vhIt.upperGlobalPos().z();
0566
0567 low_xx_global_err = vhIt.lowerGlobalPosErr().cxx();
0568 low_yy_global_err = vhIt.lowerGlobalPosErr().cyy();
0569 low_zz_global_err = vhIt.lowerGlobalPosErr().czz();
0570 low_xy_global_err = vhIt.lowerGlobalPosErr().cyx();
0571 low_zx_global_err = vhIt.lowerGlobalPosErr().czx();
0572 low_zy_global_err = vhIt.lowerGlobalPosErr().czy();
0573
0574 upp_xx_global_err = vhIt.upperGlobalPosErr().cxx();
0575 upp_yy_global_err = vhIt.upperGlobalPosErr().cyy();
0576 upp_zz_global_err = vhIt.upperGlobalPosErr().czz();
0577 upp_xy_global_err = vhIt.upperGlobalPosErr().cyx();
0578 upp_zx_global_err = vhIt.upperGlobalPosErr().czx();
0579 upp_zy_global_err = vhIt.upperGlobalPosErr().czy();
0580
0581 LogDebug("VectorHitsBuilderValidation") << "print Clusters into the VH:" << std::endl;
0582 printCluster(geomDetLower, vhIt.lowerClusterRef());
0583 LogTrace("VectorHitsBuilderValidation") << "\t global pos lower " << vhIt.lowerGlobalPos() << std::endl;
0584 LogTrace("VectorHitsBuilderValidation")
0585 << "\t global posErr lower " << vhIt.lowerGlobalPosErr().cxx() << std::endl;
0586 const GeomDetUnit* geomDetUpper = stackDet->upperDet();
0587 printCluster(geomDetUpper, vhIt.upperClusterRef());
0588 LogTrace("VectorHitsBuilderValidation") << "\t global pos upper " << vhIt.upperGlobalPos() << std::endl;
0589
0590
0591 LogDebug("VectorHitsBuilderValidation") << "comparison Clusters with sim hits ... " << std::endl;
0592 std::vector<unsigned int> clusterSimTrackIds;
0593 std::vector<unsigned int> clusterSimTrackIdsUpp;
0594 std::set<std::pair<uint32_t, EncodedEventId> > simTkIds;
0595 const GeomDetUnit* geomDetUnit_low(tkGeom_->idToDetUnit(lowerDetId));
0596 LogTrace("VectorHitsBuilderValidation") << " lowerDetID : " << lowerDetId.rawId();
0597 const GeomDetUnit* geomDetUnit_upp(tkGeom_->idToDetUnit(upperDetId));
0598 LogTrace("VectorHitsBuilderValidation") << " upperDetID : " << upperDetId.rawId();
0599
0600 for (unsigned int istr(0); istr < (*(vhIt.lowerClusterRef().cluster_phase2OT())).size(); ++istr) {
0601 uint32_t channel =
0602 Phase2TrackerDigi::pixelToChannel((*(vhIt.lowerClusterRef().cluster_phase2OT())).firstRow() + istr,
0603 (*(vhIt.lowerClusterRef().cluster_phase2OT())).column());
0604 unsigned int LowerSimTrackId(getSimTrackId(siphase2SimLinks, lowerDetId, channel));
0605 std::vector<std::pair<uint32_t, EncodedEventId> > trkid(
0606 getSimTrackIds(siphase2SimLinks, lowerDetId, channel));
0607 if (trkid.size() == 0)
0608 continue;
0609 clusterSimTrackIds.push_back(LowerSimTrackId);
0610 simTkIds.insert(trkid.begin(), trkid.end());
0611 LogTrace("VectorHitsBuilderValidation") << "LowerSimTrackId " << LowerSimTrackId << std::endl;
0612 }
0613
0614 for (const auto& iset : simTkIds) {
0615 auto ipos = mapping.find(iset);
0616 if (ipos != mapping.end()) {
0617 LogTrace("VectorHitsBuilderValidation") << "lower cluster in detid: " << lowerDetId.rawId()
0618 << " from tp: " << ipos->second.key() << " " << iset.first;
0619 LogTrace("VectorHitsBuilderValidation") << "with pt(): " << (*ipos->second).pt();
0620 low_tp_id = ipos->second.key();
0621 vh_sim_trackPt = (*ipos->second).pt();
0622 }
0623 }
0624
0625 simTkIds.clear();
0626 for (unsigned int istr(0); istr < (*(vhIt.upperClusterRef().cluster_phase2OT())).size(); ++istr) {
0627 uint32_t channel =
0628 Phase2TrackerDigi::pixelToChannel((*(vhIt.upperClusterRef().cluster_phase2OT())).firstRow() + istr,
0629 (*(vhIt.upperClusterRef().cluster_phase2OT())).column());
0630 unsigned int UpperSimTrackId(getSimTrackId(siphase2SimLinks, upperDetId, channel));
0631 std::vector<std::pair<uint32_t, EncodedEventId> > trkid(
0632 getSimTrackIds(siphase2SimLinks, upperDetId, channel));
0633 if (trkid.size() == 0)
0634 continue;
0635 clusterSimTrackIdsUpp.push_back(UpperSimTrackId);
0636 simTkIds.insert(trkid.begin(), trkid.end());
0637 LogTrace("VectorHitsBuilderValidation") << "UpperSimTrackId " << UpperSimTrackId << std::endl;
0638 }
0639
0640 for (const auto& iset : simTkIds) {
0641 auto ipos = mapping.find(iset);
0642 if (ipos != mapping.end()) {
0643 LogTrace("VectorHitsBuilderValidation")
0644 << "upper cluster in detid: " << upperDetId.rawId() << " from tp: " << ipos->second.key() << " "
0645 << iset.first << std::endl;
0646 upp_tp_id = ipos->second.key();
0647 }
0648 }
0649
0650 std::pair<bool, uint32_t> istrue = isTrue(vhIt, siphase2SimLinks, detId);
0651 vh_isTrue = 0;
0652 if (istrue.first) {
0653 vh_isTrue = 1;
0654 }
0655
0656
0657 unsigned int totalSimHits(0);
0658 unsigned int primarySimHits(0);
0659 unsigned int otherSimHits(0);
0660
0661 for (const auto& hitIt : *simHitsRaw) {
0662 if (hitIt.detUnitId() == geomDetLower->geographicalId()) {
0663
0664 if (std::find(clusterSimTrackIds.begin(), clusterSimTrackIds.end(), hitIt.trackId()) !=
0665 clusterSimTrackIds.end()) {
0666 Local3DPoint localPosHit(hitIt.localPosition());
0667 sim_x_local = localPosHit.x();
0668 sim_y_local = localPosHit.y();
0669
0670 deltaXVHSimHits = vh_x_local - sim_x_local;
0671 deltaYVHSimHits = vh_y_local - sim_y_local;
0672
0673 Global3DPoint globalPosHit = geomDetLower->surface().toGlobal(localPosHit);
0674 sim_x_global = globalPosHit.x();
0675 sim_y_global = globalPosHit.y();
0676 sim_z_global = globalPosHit.z();
0677
0678 histogramLayer->second.deltaXVHSimHits[0]->Fill(vh_x_local - sim_x_local);
0679 histogramLayer->second.deltaYVHSimHits[0]->Fill(vh_y_local - sim_y_local);
0680
0681
0682 if (layer == 1 || layer == 2 || layer == 3) {
0683 histogramLayer->second.deltaXVHSimHits[1]->Fill(vh_x_local - sim_x_local);
0684 histogramLayer->second.deltaYVHSimHits[1]->Fill(vh_y_local - sim_y_local);
0685 }
0686
0687 else if (layer == 4 || layer == 5 || layer == 6) {
0688 histogramLayer->second.deltaXVHSimHits[2]->Fill(vh_x_local - sim_x_local);
0689 histogramLayer->second.deltaYVHSimHits[2]->Fill(vh_y_local - sim_y_local);
0690 }
0691
0692 ++totalSimHits;
0693
0694 std::map<unsigned int, SimTrack>::const_iterator simTrackIt(simTracks.find(hitIt.trackId()));
0695 if (simTrackIt == simTracks.end())
0696 continue;
0697
0698
0699 processType = hitIt.processType();
0700
0701 if (simTrackIt->second.vertIndex() == 0 and
0702 (processType == 2 || processType == 7 || processType == 9 || processType == 11 || processType == 13 ||
0703 processType == 15)) {
0704 histogramLayer->second.deltaXVHSimHits_P[0]->Fill(vh_x_local - sim_x_local);
0705 histogramLayer->second.deltaYVHSimHits_P[0]->Fill(vh_y_local - sim_y_local);
0706
0707
0708 if (layer == 1 || layer == 2 || layer == 3) {
0709 histogramLayer->second.deltaXVHSimHits_P[1]->Fill(vh_x_local - sim_x_local);
0710 histogramLayer->second.deltaYVHSimHits_P[1]->Fill(vh_y_local - sim_y_local);
0711 }
0712
0713 else if (layer == 4 || layer == 5 || layer == 6) {
0714 histogramLayer->second.deltaXVHSimHits_P[2]->Fill(vh_x_local - sim_x_local);
0715 histogramLayer->second.deltaYVHSimHits_P[2]->Fill(vh_y_local - sim_y_local);
0716 }
0717
0718 ++primarySimHits;
0719 }
0720
0721 otherSimHits = totalSimHits - primarySimHits;
0722
0723 histogramLayer->second.totalSimHits->Fill(totalSimHits);
0724 histogramLayer->second.primarySimHits->Fill(primarySimHits);
0725 histogramLayer->second.otherSimHits->Fill(otherSimHits);
0726 }
0727 }
0728 }
0729
0730 nVHsTot++;
0731
0732
0733
0734
0735 multiplicity = DSViter.size();
0736 if (DSViter.size() > 1) {
0737 LogTrace("VectorHitsBuilderValidation") << " not filling if more than 1 VH has been produced";
0738 width = -100;
0739 deltaXlocal = -100;
0740 tree_->Fill();
0741 continue;
0742 }
0743
0744
0745 GlobalPoint center(0.0, 0.0, 0.0);
0746 curvature = vh.curvature();
0747 phi = vh.phi();
0748 QOverPT = vh.transverseMomentum(magField_->inTesla(center).z());
0749 QOverP = vh.momentum(magField_->inTesla(center).z());
0750 histogramLayer->second.curvature->Fill(curvature);
0751
0752
0753
0754 auto&& lparamsUpp = cpe_->localParameters(*vhIt.upperClusterRef().cluster_phase2OT(), *geomDetUnit_upp);
0755 LogTrace("VectorHitsBuilderValidation") << " upper local pos (in its system of reference):" << lparamsUpp.first;
0756 Global3DPoint gparamsUpp = geomDetUnit_upp->surface().toGlobal(lparamsUpp.first);
0757 LogTrace("VectorHitsBuilderValidation") << " upper global pos :" << gparamsUpp;
0758 Local3DPoint lparamsUppInLow = geomDetUnit_low->surface().toLocal(gparamsUpp);
0759 LogTrace("VectorHitsBuilderValidation") << " upper local pos (in low system of reference):" << lparamsUppInLow;
0760 auto&& lparamsLow = cpe_->localParameters(*vhIt.lowerClusterRef().cluster_phase2OT(), *geomDetUnit_low);
0761 LogTrace("VectorHitsBuilderValidation") << " lower local pos (in its system of reference):" << lparamsLow.first;
0762 Global3DPoint gparamsLow = geomDetUnit_low->surface().toGlobal(lparamsLow.first);
0763 LogTrace("VectorHitsBuilderValidation") << " lower global pos :" << gparamsLow;
0764
0765 deltaXlocal = lparamsUppInLow.x() - lparamsLow.first.x();
0766 histogramLayer->second.deltaXlocal->Fill(deltaXlocal);
0767 LogTrace("VectorHitsBuilderValidation") << " deltaXlocal : " << deltaXlocal;
0768
0769 double parallCorr = 0.0;
0770
0771 Global3DPoint origin(0, 0, 0);
0772 GlobalVector gV = gparamsLow - origin;
0773 LocalVector lV = geomDetUnit_low->surface().toLocal(gV);
0774 LocalVector lV_norm = lV / lV.z();
0775 parallCorr = lV_norm.x() * lparamsUppInLow.z();
0776 LogTrace("VectorHitsBuilderValidation") << " parallalex correction:" << parallCorr;
0777
0778 double lpos_upp_corr = 0.0;
0779 double lpos_low_corr = 0.0;
0780 if (lparamsUpp.first.x() > lparamsLow.first.x()) {
0781 if (lparamsUpp.first.x() > 0) {
0782 lpos_low_corr = lparamsLow.first.x();
0783 lpos_upp_corr = lparamsUpp.first.x() - fabs(parallCorr);
0784 }
0785 if (lparamsUpp.first.x() < 0) {
0786 lpos_low_corr = lparamsLow.first.x() + fabs(parallCorr);
0787 lpos_upp_corr = lparamsUpp.first.x();
0788 }
0789 } else if (lparamsUpp.first.x() < lparamsLow.first.x()) {
0790 if (lparamsUpp.first.x() > 0) {
0791 lpos_low_corr = lparamsLow.first.x() - fabs(parallCorr);
0792 lpos_upp_corr = lparamsUpp.first.x();
0793 }
0794 if (lparamsUpp.first.x() < 0) {
0795 lpos_low_corr = lparamsLow.first.x();
0796 lpos_upp_corr = lparamsUpp.first.x() + fabs(parallCorr);
0797 }
0798 } else {
0799 if (lparamsUpp.first.x() > 0) {
0800 lpos_upp_corr = lparamsUpp.first.x() - fabs(parallCorr);
0801 lpos_low_corr = lparamsLow.first.x();
0802 }
0803 if (lparamsUpp.first.x() < 0) {
0804 lpos_upp_corr = lparamsUpp.first.x() + fabs(parallCorr);
0805 lpos_low_corr = lparamsLow.first.x();
0806 }
0807 }
0808
0809 LogDebug("VectorHitsBuilderValidation") << " \t local pos upper corrected (x):" << lpos_upp_corr << std::endl;
0810 LogDebug("VectorHitsBuilderValidation") << " \t local pos lower corrected (x):" << lpos_low_corr << std::endl;
0811
0812 width = lpos_low_corr - lpos_upp_corr;
0813 histogramLayer->second.width->Fill(width);
0814 LogTrace("VectorHitsBuilderValidation") << " width:" << width;
0815
0816 tree_->Fill();
0817
0818 }
0819
0820 }
0821
0822 if (nVHsPS)
0823 histogramLayer->second.numberVHsPS->Fill(nVHsPS);
0824 if (nVHs2S)
0825 histogramLayer->second.numberVHs2S->Fill(nVHs2S);
0826 LogTrace("VectorHitsBuilderValidation")
0827 << "nVHsPS for this layer : " << nVHsPS << ", nVHs2S for this layer : " << nVHs2S << std::endl;
0828 }
0829
0830 CreateVHsXYGraph(glVHs, dirVHs);
0831 CreateVHsRZGraph(glVHs, dirVHs);
0832
0833 int VHaccTrue = 0.0;
0834 int VHaccFalse = 0.0;
0835 int VHrejTrue = 0.0;
0836 int VHrejFalse = 0.0;
0837 int VHaccTrue_signal = 0.0;
0838 int VHrejTrue_signal = 0.0;
0839
0840
0841 for (const auto& DSViter : *vhsAcc) {
0842 unsigned int rawid(DSViter.detId());
0843 DetId detId(rawid);
0844 int layerAcc = getLayerNumber(detId);
0845 LogTrace("VectorHitsBuilderValidation") << "acc Layer: " << layerAcc << " det id" << rawid << std::endl;
0846 for (const auto& vhIt : DSViter) {
0847 if (vhIt.isValid()) {
0848 VHaccLayer_->Fill(layerAcc);
0849 VHacc++;
0850
0851
0852 std::pair<bool, uint32_t> istrue = isTrue(vhIt, siphase2SimLinks, detId);
0853 if (istrue.first) {
0854 LogTrace("VectorHitsBuilderValidation") << "this vectorhit is a 'true' vhit.";
0855 VHaccTrueLayer_->Fill(layerAcc);
0856 VHaccTrue++;
0857
0858
0859 std::map<unsigned int, SimTrack>::const_iterator simTrackIt(simTracks.find(istrue.second));
0860 if (simTrackIt == simTracks.end())
0861 continue;
0862 LogTrace("VectorHitsBuilderValidation") << "this vectorhit is associated with SimTrackId: " << istrue.second;
0863 LogTrace("VectorHitsBuilderValidation") << "the SimTrack has pt: " << simTrackIt->second.momentum().pt();
0864 if (simTrackIt->second.momentum().pt() > 1) {
0865 VHaccTrue_signal_Layer_->Fill(layerAcc);
0866 LogTrace("VectorHitsBuilderValidation") << "the vectorhit belongs to signal";
0867 VHaccTrue_signal++;
0868 }
0869
0870 } else {
0871 LogTrace("VectorHitsBuilderValidation") << "this vectorhit is a 'false' vhit.";
0872 VHaccFalse++;
0873 }
0874 }
0875 }
0876 }
0877
0878 for (const auto& DSViter : *vhsRej) {
0879 unsigned int rawid(DSViter.detId());
0880 DetId detId(rawid);
0881 int layerRej = getLayerNumber(detId);
0882 LogTrace("VectorHitsBuilderValidation") << "rej Layer: " << layerRej << " det id" << rawid << std::endl;
0883 for (const auto& vhIt : DSViter) {
0884 VHrejLayer_->Fill(layerRej);
0885 VHrej++;
0886
0887
0888 std::pair<bool, uint32_t> istrue = isTrue(vhIt, siphase2SimLinks, detId);
0889 if (istrue.first) {
0890 LogTrace("VectorHitsBuilderValidation") << "this vectorhit is a 'true' vhit.";
0891 VHrejTrueLayer_->Fill(layerRej);
0892 VHrejTrue++;
0893
0894
0895 std::map<unsigned int, SimTrack>::const_iterator simTrackIt(simTracks.find(istrue.second));
0896 if (simTrackIt == simTracks.end())
0897 continue;
0898 LogTrace("VectorHitsBuilderValidation") << "this vectorhit is associated with SimTrackId: " << istrue.second;
0899 LogTrace("VectorHitsBuilderValidation") << "the SimTrack has pt: " << simTrackIt->second.momentum().pt();
0900 if (simTrackIt->second.momentum().pt() > 1) {
0901 VHrejTrue_signal_Layer_->Fill(layerRej);
0902 LogTrace("VectorHitsBuilderValidation") << "the vectorhit belongs to signal";
0903 VHrejTrue_signal++;
0904 }
0905
0906 } else {
0907 LogTrace("VectorHitsBuilderValidation") << "this vectorhit is a 'false' vhit.";
0908 VHrejFalse++;
0909 }
0910 }
0911 }
0912
0913 int VHtot = VHacc + VHrej;
0914 LogTrace("VectorHitsBuilderValidation")
0915 << "VH total: " << VHtot << " with " << VHacc << " VHs accepted and " << VHrej << " VHs rejected.";
0916 LogTrace("VectorHitsBuilderValidation")
0917 << "of the VH accepted, there are " << VHaccTrue << " true and " << VHaccFalse << " false.";
0918 LogTrace("VectorHitsBuilderValidation")
0919 << "of the VH rejected, there are " << VHrejTrue << " true and " << VHrejFalse << " false.";
0920 LogTrace("VectorHitsBuilderValidation")
0921 << "of the true VH , there are " << VHaccTrue_signal << " accepted belonging to signal and "
0922 << VHrejTrue_signal << " rejected belonging to signal.";
0923
0924
0925 }
0926
0927
0928 std::pair<bool, uint32_t> VectorHitsBuilderValidation::isTrue(
0929 const VectorHit vh, const edm::Handle<edm::DetSetVector<PixelDigiSimLink> >& siphase2SimLinks, DetId& detId) const {
0930 const GeomDet* geomDet(tkGeom_->idToDet(detId));
0931 const StackGeomDet* stackDet = dynamic_cast<const StackGeomDet*>(geomDet);
0932 const GeomDetUnit* geomDetLower = stackDet->lowerDet();
0933 const GeomDetUnit* geomDetUpper = stackDet->upperDet();
0934
0935 std::vector<unsigned int> lowClusterSimTrackIds;
0936
0937 for (unsigned int istr(0); istr < (*(vh.lowerClusterRef().cluster_phase2OT())).size(); ++istr) {
0938 uint32_t channel = Phase2TrackerDigi::pixelToChannel((*(vh.lowerClusterRef().cluster_phase2OT())).firstRow() + istr,
0939 (*(vh.lowerClusterRef().cluster_phase2OT())).column());
0940 DetId detIdCluster = geomDetLower->geographicalId();
0941 unsigned int simTrackId(getSimTrackId(siphase2SimLinks, detIdCluster, channel));
0942 LogTrace("VectorHitsBuilderValidation") << "LowerSimTrackId " << simTrackId << std::endl;
0943 std::vector<std::pair<uint32_t, EncodedEventId> > trkid(getSimTrackIds(siphase2SimLinks, detIdCluster, channel));
0944 if (trkid.size() == 0)
0945 continue;
0946 lowClusterSimTrackIds.push_back(simTrackId);
0947 }
0948
0949 std::vector<unsigned int>::iterator it_simTrackUpper;
0950
0951 for (unsigned int istr(0); istr < (*(vh.upperClusterRef().cluster_phase2OT())).size(); ++istr) {
0952 uint32_t channel = Phase2TrackerDigi::pixelToChannel((*(vh.upperClusterRef().cluster_phase2OT())).firstRow() + istr,
0953 (*(vh.upperClusterRef().cluster_phase2OT())).column());
0954 DetId detIdCluster = geomDetUpper->geographicalId();
0955 unsigned int simTrackId(getSimTrackId(siphase2SimLinks, detIdCluster, channel));
0956 LogTrace("VectorHitsBuilderValidation") << "UpperSimTrackId " << simTrackId << std::endl;
0957 std::vector<std::pair<uint32_t, EncodedEventId> > trkid(getSimTrackIds(siphase2SimLinks, detIdCluster, channel));
0958 if (trkid.size() == 0)
0959 continue;
0960 it_simTrackUpper = std::find(lowClusterSimTrackIds.begin(), lowClusterSimTrackIds.end(), simTrackId);
0961 if (it_simTrackUpper != lowClusterSimTrackIds.end()) {
0962 LogTrace("VectorHitsBuilderValidation") << " UpperSimTrackId found in lowClusterSimTrackIds ";
0963 return std::make_pair(true, simTrackId);
0964 }
0965 }
0966 return std::make_pair(false, 0);
0967 }
0968
0969
0970 std::map<unsigned int, VHHistos>::iterator VectorHitsBuilderValidation::createLayerHistograms(unsigned int ival) {
0971 std::ostringstream fname1, fname2;
0972
0973 edm::Service<TFileService> fs;
0974 fs->file().cd("/");
0975
0976 std::string tag;
0977 unsigned int id;
0978 if (ival < 100) {
0979 id = ival;
0980 fname1 << "Barrel";
0981 fname2 << "Layer_" << id;
0982 tag = "_layer_";
0983 } else {
0984 int side = ival / 100;
0985 id = ival - side * 100;
0986 fname1 << "EndCap_Side_" << side;
0987 fname2 << "Disc_" << id;
0988 tag = "_disc_";
0989 }
0990
0991 TFileDirectory td1 = fs->mkdir(fname1.str().c_str());
0992 TFileDirectory td = td1.mkdir(fname2.str().c_str());
0993
0994 VHHistos local_histos;
0995
0996 std::ostringstream histoName;
0997
0998
0999
1000
1001
1002 histoName.str("");
1003 histoName << "Number_VHs_PS" << tag.c_str() << id;
1004 local_histos.numberVHsPS = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 20, 0., 20.);
1005 local_histos.numberVHsPS->SetFillColor(kAzure + 7);
1006
1007 histoName.str("");
1008 histoName << "Number_VHs_2S" << tag.c_str() << id;
1009 local_histos.numberVHs2S = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 20, 0., 20.);
1010 local_histos.numberVHs2S->SetFillColor(kOrange - 3);
1011
1012 histoName.str("");
1013 histoName << "Number_VHs_Mixed" << tag.c_str() << id;
1014 local_histos.numberVHsMixed = td.make<THStack>(histoName.str().c_str(), histoName.str().c_str());
1015 local_histos.numberVHsMixed->Add(local_histos.numberVHsPS);
1016 local_histos.numberVHsMixed->Add(local_histos.numberVHs2S);
1017
1018
1019
1020
1021
1022 histoName.str("");
1023 histoName << "Local_Position_XY_Mixed" << tag.c_str() << id;
1024 local_histos.localPosXY[0] = td.make<TGraph>();
1025 local_histos.localPosXY[0]->SetName(histoName.str().c_str());
1026
1027 histoName.str("");
1028 histoName << "Local_Position_XY_PS" << tag.c_str() << id;
1029 local_histos.localPosXY[1] = td.make<TGraph>();
1030 local_histos.localPosXY[1]->SetName(histoName.str().c_str());
1031
1032 histoName.str("");
1033 histoName << "Local_Position_XY_2S" << tag.c_str() << id;
1034 local_histos.localPosXY[2] = td.make<TGraph>();
1035 local_histos.localPosXY[2]->SetName(histoName.str().c_str());
1036
1037 histoName.str("");
1038 histoName << "Global_Position_XY_Mixed" << tag.c_str() << id;
1039 local_histos.globalPosXY[0] = td.make<TGraph>();
1040 local_histos.globalPosXY[0]->SetName(histoName.str().c_str());
1041
1042 histoName.str("");
1043 histoName << "Global_Position_XY_PS" << tag.c_str() << id;
1044 local_histos.globalPosXY[1] = td.make<TGraph>();
1045 local_histos.globalPosXY[1]->SetName(histoName.str().c_str());
1046
1047 histoName.str("");
1048 histoName << "Global_Position_XY_2S" << tag.c_str() << id;
1049 local_histos.globalPosXY[2] = td.make<TGraph>();
1050 local_histos.globalPosXY[2]->SetName(histoName.str().c_str());
1051
1052
1053
1054
1055
1056 histoName.str("");
1057 histoName << "Delta_X_VH_SimHits_Mixed" << tag.c_str() << id;
1058 local_histos.deltaXVHSimHits[0] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
1059
1060 histoName.str("");
1061 histoName << "Delta_X_VH_SimHits_PS" << tag.c_str() << id;
1062 local_histos.deltaXVHSimHits[1] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
1063
1064 histoName.str("");
1065 histoName << "Delta_X_VH_SimHits_2S" << tag.c_str() << id;
1066 local_histos.deltaXVHSimHits[2] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
1067
1068 histoName.str("");
1069 histoName << "Delta_Y_VH_SimHits_Mixed" << tag.c_str() << id;
1070 local_histos.deltaYVHSimHits[0] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
1071
1072 histoName.str("");
1073 histoName << "Delta_Y_VH_SimHits_PS" << tag.c_str() << id;
1074 local_histos.deltaYVHSimHits[1] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
1075
1076 histoName.str("");
1077 histoName << "Delta_Y_VH_SimHits_2S" << tag.c_str() << id;
1078 local_histos.deltaYVHSimHits[2] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
1079
1080
1081
1082
1083
1084 histoName.str("");
1085 histoName << "Delta_X_VH_SimHits_Mixed_P" << tag.c_str() << id;
1086 local_histos.deltaXVHSimHits_P[0] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
1087
1088 histoName.str("");
1089 histoName << "Delta_X_VH_SimHits_PS_P" << tag.c_str() << id;
1090 local_histos.deltaXVHSimHits_P[1] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
1091
1092 histoName.str("");
1093 histoName << "Delta_X_VH_SimHits_2S_P" << tag.c_str() << id;
1094 local_histos.deltaXVHSimHits_P[2] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
1095
1096 histoName.str("");
1097 histoName << "Delta_Y_VH_SimHits_Mixed_P" << tag.c_str() << id;
1098 local_histos.deltaYVHSimHits_P[0] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
1099
1100 histoName.str("");
1101 histoName << "Delta_Y_VH_SimHits_PS_P" << tag.c_str() << id;
1102 local_histos.deltaYVHSimHits_P[1] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
1103
1104 histoName.str("");
1105 histoName << "Delta_Y_VH_SimHits_2S_P" << tag.c_str() << id;
1106 local_histos.deltaYVHSimHits_P[2] = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 100, 0., 0.);
1107
1108
1109
1110
1111
1112 histoName.str("");
1113 histoName << "Total_Digis" << tag.c_str() << id;
1114 local_histos.totalSimHits = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 10, 0., 10.);
1115
1116 histoName.str("");
1117 histoName << "Primary_Digis" << tag.c_str() << id;
1118 local_histos.primarySimHits = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 10, 0., 10.);
1119
1120 histoName.str("");
1121 histoName << "Other_Digis" << tag.c_str() << id;
1122 local_histos.otherSimHits = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 10, 0., 10.);
1123
1124
1125
1126
1127
1128 histoName.str("");
1129 histoName << "DeltaXlocal_clusters" << tag.c_str() << id;
1130 local_histos.deltaXlocal = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 200, -0.4, 0.4);
1131 histoName.str("");
1132 histoName << "Width" << tag.c_str() << id;
1133 local_histos.width = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 200, -0.4, 0.4);
1134 histoName.str("");
1135 histoName << "Curvature" << tag.c_str() << id;
1136 local_histos.curvature = td.make<TH1F>(histoName.str().c_str(), histoName.str().c_str(), 200, -0.4, 0.4);
1137
1138 std::pair<std::map<unsigned int, VHHistos>::iterator, bool> insertedIt(
1139 histograms_.insert(std::make_pair(ival, local_histos)));
1140 fs->file().cd("/");
1141
1142 return insertedIt.first;
1143 }
1144
1145 void VectorHitsBuilderValidation::CreateVHsXYGraph(const std::vector<Global3DPoint> glVHs,
1146 const std::vector<Global3DVector> dirVHs) {
1147 if (glVHs.size() != dirVHs.size()) {
1148 std::cout << "Cannot fullfil the graphs for this event. Return." << std::endl;
1149 return;
1150 }
1151
1152
1153
1154 for (unsigned int nVH = 0; nVH < glVHs.size(); nVH++) {
1155
1156 if ((fabs(dirVHs.at(nVH).x()) < 10e-5) && (fabs(dirVHs.at(nVH).y()) < 10e-5)) {
1157 continue;
1158
1159 } else {
1160 }
1161 }
1162
1163 return;
1164 }
1165
1166 void VectorHitsBuilderValidation::CreateVHsRZGraph(const std::vector<Global3DPoint> glVHs,
1167 const std::vector<Global3DVector> dirVHs) {
1168 if (glVHs.size() != dirVHs.size()) {
1169 std::cout << "Cannot fullfil the graphs for this event. Return." << std::endl;
1170 return;
1171 }
1172
1173 return;
1174 }
1175
1176 void VectorHitsBuilderValidation::CreateWindowCorrGraph() {
1177
1178
1179
1180 if (!tkGeom_->isThere(GeomDetEnumerators::P2OTB) && !tkGeom_->isThere(GeomDetEnumerators::P2OTEC))
1181 return;
1182
1183 for (auto det : tkGeom_->detsTOB()) {
1184 ParallaxCorrectionRZ_->Fill(det->position().z(), det->position().perp(), 5.);
1185 }
1186 for (auto det : tkGeom_->detsTID()) {
1187 ParallaxCorrectionRZ_->Fill(det->position().z(), det->position().perp(), 10.);
1188 }
1189 ParallaxCorrectionRZ_->Fill(0., 0., 5.);
1190 return;
1191 }
1192
1193 unsigned int VectorHitsBuilderValidation::getLayerNumber(const DetId& detid) {
1194 if (detid.det() == DetId::Tracker) {
1195 if (detid.subdetId() == StripSubdetector::TOB)
1196 return (tkTopo_->layer(detid));
1197 else if (detid.subdetId() == StripSubdetector::TID)
1198 return (100 * tkTopo_->side(detid) + tkTopo_->layer(detid));
1199 else
1200 return 999;
1201 }
1202 return 999;
1203 }
1204
1205 unsigned int VectorHitsBuilderValidation::getModuleNumber(const DetId& detid) { return (tkTopo_->module(detid)); }
1206
1207 std::vector<std::pair<uint32_t, EncodedEventId> > VectorHitsBuilderValidation::getSimTrackIds(
1208 const edm::Handle<edm::DetSetVector<PixelDigiSimLink> >& simLinks, const DetId& detId, uint32_t channel) const {
1209 std::vector<std::pair<uint32_t, EncodedEventId> > simTrkId;
1210 auto isearch = simLinks->find(detId);
1211 if (isearch != simLinks->end()) {
1212
1213 edm::DetSet<PixelDigiSimLink> link_detset = (*isearch);
1214 for (const auto& it : link_detset.data) {
1215 if (channel == it.channel())
1216 simTrkId.push_back(std::make_pair(it.SimTrackId(), it.eventId()));
1217 }
1218 }
1219 return simTrkId;
1220 }
1221
1222 unsigned int VectorHitsBuilderValidation::getSimTrackId(
1223 const edm::Handle<edm::DetSetVector<PixelDigiSimLink> >& pixelSimLinks,
1224 const DetId& detId,
1225 unsigned int channel) const {
1226 edm::DetSetVector<PixelDigiSimLink>::const_iterator DSViter(pixelSimLinks->find(detId));
1227 if (DSViter == pixelSimLinks->end())
1228 return 0;
1229 for (const auto& it : DSViter->data) {
1230 if (channel == it.channel())
1231 return it.SimTrackId();
1232 }
1233 return 0;
1234 }
1235
1236 void VectorHitsBuilderValidation::printCluster(const GeomDetUnit* geomDetUnit, const OmniClusterRef cluster) {
1237 if (!geomDetUnit)
1238 return;
1239
1240 const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>(geomDetUnit);
1241 const PixelTopology& topol = theGeomDet->specificTopology();
1242
1243 unsigned int layer = getLayerNumber(geomDetUnit->geographicalId());
1244 unsigned int module = getModuleNumber(geomDetUnit->geographicalId());
1245 LogTrace("VectorHitsBuilderValidation") << "Layer:" << layer << std::endl;
1246 if (topol.ncolumns() == 32)
1247 LogTrace("VectorHitsBuilderValidation") << "Pixel cluster with detId:" << geomDetUnit->geographicalId().rawId()
1248 << "(module:" << module << ") " << std::endl;
1249 else if (topol.ncolumns() == 2)
1250 LogTrace("VectorHitsBuilderValidation") << "Strip cluster with detId:" << geomDetUnit->geographicalId().rawId()
1251 << "(module:" << module << ") " << std::endl;
1252 else
1253 std::cout << "no module?!" << std::endl;
1254 LogTrace("VectorHitsBuilderValidation")
1255 << "with pitch:" << topol.pitch().first << " , " << topol.pitch().second << std::endl;
1256 LogTrace("VectorHitsBuilderValidation") << " and width:" << theGeomDet->surface().bounds().width()
1257 << " , lenght:" << theGeomDet->surface().bounds().length() << std::endl;
1258
1259 auto&& lparams = cpe_->localParameters(*cluster.cluster_phase2OT(), *theGeomDet);
1260
1261 LogTrace("VectorHitsBuilderValidation")
1262 << "\t local pos " << lparams.first << "with err " << lparams.second << std::endl;
1263
1264 return;
1265 }
1266
1267 void VectorHitsBuilderValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1268 edm::ParameterSetDescription desc;
1269 desc.add<std::string>("src", "siPhase2Clusters");
1270 desc.add<edm::InputTag>("links", edm::InputTag("simSiPixelDigis", "Tracker"));
1271 desc.add<edm::InputTag>("VH_acc", edm::InputTag("siPhase2VectorHits", "accepted"));
1272 desc.add<edm::InputTag>("VH_rej", edm::InputTag("siPhase2VectorHits", "rejected"));
1273 desc.add<edm::ESInputTag>("CPE", edm::ESInputTag("phase2StripCPEESProducer", "Phase2StripCPE"));
1274 desc.add<edm::InputTag>("trackingParticleSrc", edm::InputTag("mix", "MergedTrackTruth"));
1275 descriptions.add("vectorHitsBuilderValidation", desc);
1276 }
1277
1278 DEFINE_FWK_MODULE(VectorHitsBuilderValidation);
