File indexing completed on 2024-04-06 12:28:04
0001 #include "TestTrackHits.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0004 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0005 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0006 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0007 #include <TDirectory.h>
0008
0009 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0010 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0011 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0012
0013 typedef TrajectoryStateOnSurface TSOS;
0014 typedef TransientTrackingRecHit::ConstRecHitPointer CTTRHp;
0015 using namespace std;
0016 using namespace edm;
0017
0018 TestTrackHits::TestTrackHits(const edm::ParameterSet& iConfig) : trackerHitAssociatorConfig_(consumesCollector()) {
0019 LogTrace("TestTrackHits") << iConfig;
0020 propagatorName = iConfig.getParameter<std::string>("Propagator");
0021 builderName = iConfig.getParameter<std::string>("TTRHBuilder");
0022 srcName = iConfig.getParameter<std::string>("src");
0023 tpName = iConfig.getParameter<std::string>("tpname");
0024 updatorName = iConfig.getParameter<std::string>("updator");
0025 out = iConfig.getParameter<std::string>("out");
0026 theGToken = esConsumes<edm::Transition::BeginRun>();
0027 theMFToken = esConsumes<edm::Transition::BeginRun>();
0028 thePropagatorToken = esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", propagatorName));
0029 theBuilderToken = esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", builderName));
0030 theUpdatorToken = esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", updatorName));
0031 theTopoToken = esConsumes();
0032
0033
0034
0035
0036
0037
0038
0039
0040 }
0041
0042 TestTrackHits::~TestTrackHits() {}
0043
0044 void TestTrackHits::beginRun(edm::Run const& run, const edm::EventSetup& iSetup) {
0045 theG = iSetup.getHandle(theGToken);
0046 theMF = iSetup.getHandle(theMFToken);
0047 thePropagator = iSetup.getHandle(thePropagatorToken);
0048 theBuilder = iSetup.getHandle(theBuilderToken);
0049 theUpdator = iSetup.getHandle(theUpdatorToken);
0050
0051 file = new TFile(out.c_str(), "recreate");
0052 for (int i = 0; i != 6; i++)
0053 for (int j = 0; j != 9; j++) {
0054 if (i == 0 && j > 2)
0055 break;
0056 if (i == 1 && j > 1)
0057 break;
0058 if (i == 2 && j > 3)
0059 break;
0060 if (i == 3 && j > 2)
0061 break;
0062 if (i == 4 && j > 5)
0063 break;
0064 if (i == 5 && j > 8)
0065 break;
0066
0067 title.str("");
0068 title << "Chi2Increment_" << i + 1 << "-" << j + 1;
0069 hChi2Increment[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, 0, 100);
0070 title.str("");
0071 title << "Chi2IncrementVsEta_" << i + 1 << "-" << j + 1;
0072 hChi2IncrementVsEta[title.str()] =
0073 new TH2F(title.str().c_str(), title.str().c_str(), 50, -2.5, 2.5, 1000, 0, 100);
0074 title.str("");
0075 title << "Chi2GoodHit_" << i + 1 << "-" << j + 1;
0076 hChi2GoodHit[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, 0, 100);
0077 title.str("");
0078 title << "Chi2BadHit_" << i + 1 << "-" << j + 1;
0079 hChi2BadHit[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, 0, 100);
0080 title.str("");
0081 title << "Chi2DeltaHit_" << i + 1 << "-" << j + 1;
0082 hChi2DeltaHit[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, 0, 100);
0083 title.str("");
0084 title << "Chi2NSharedHit_" << i + 1 << "-" << j + 1;
0085 hChi2NSharedHit[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, 0, 100);
0086 title.str("");
0087 title << "Chi2SharedHit_" << i + 1 << "-" << j + 1;
0088 hChi2SharedHit[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, 0, 100);
0089
0090 title.str("");
0091 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_ts";
0092 hPullGP_X_ts[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0093 title.str("");
0094 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_ts";
0095 hPullGP_Y_ts[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0096 title.str("");
0097 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_ts";
0098 hPullGP_Z_ts[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0099
0100 title.str("");
0101 title << "PullGM_X_" << i + 1 << "-" << j + 1 << "_ts";
0102 hPullGM_X_ts[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0103 title.str("");
0104 title << "PullGM_Y_" << i + 1 << "-" << j + 1 << "_ts";
0105 hPullGM_Y_ts[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0106 title.str("");
0107 title << "PullGM_Z_" << i + 1 << "-" << j + 1 << "_ts";
0108 hPullGM_Z_ts[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0109
0110 title.str("");
0111 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_tr";
0112 hPullGP_X_tr[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0113 title.str("");
0114 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_tr";
0115 hPullGP_Y_tr[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0116 title.str("");
0117 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_tr";
0118 hPullGP_Z_tr[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0119
0120 title.str("");
0121 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_rs";
0122 hPullGP_X_rs[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0123 title.str("");
0124 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_rs";
0125 hPullGP_Y_rs[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0126 title.str("");
0127 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_rs";
0128 hPullGP_Z_rs[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0129
0130 if (((i == 2 || i == 4) && (j == 0 || j == 1)) || (i == 3 || i == 5)) {
0131
0132 title.str("");
0133 title << "Chi2Increment_mono_" << i + 1 << "-" << j + 1;
0134 hChi2Increment_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, 0, 100);
0135
0136 title.str("");
0137 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_ts_mono";
0138 hPullGP_X_ts_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0139 title.str("");
0140 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_ts_mono";
0141 hPullGP_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0142 title.str("");
0143 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_ts_mono";
0144 hPullGP_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0145
0146 title.str("");
0147 title << "PullGM_X_" << i + 1 << "-" << j + 1 << "_ts_mono";
0148 hPullGM_X_ts_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0149 title.str("");
0150 title << "PullGM_Y_" << i + 1 << "-" << j + 1 << "_ts_mono";
0151 hPullGM_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0152 title.str("");
0153 title << "PullGM_Z_" << i + 1 << "-" << j + 1 << "_ts_mono";
0154 hPullGM_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0155
0156 title.str("");
0157 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_tr_mono";
0158 hPullGP_X_tr_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0159 title.str("");
0160 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_tr_mono";
0161 hPullGP_Y_tr_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0162 title.str("");
0163 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_tr_mono";
0164 hPullGP_Z_tr_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0165
0166 title.str("");
0167 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_rs_mono";
0168 hPullGP_X_rs_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0169 title.str("");
0170 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_rs_mono";
0171 hPullGP_Y_rs_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0172 title.str("");
0173 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_rs_mono";
0174 hPullGP_Z_rs_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0175
0176
0177 title.str("");
0178 title << "Chi2Increment_stereo_" << i + 1 << "-" << j + 1;
0179 hChi2Increment_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, 0, 100);
0180
0181 title.str("");
0182 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0183 hPullGP_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0184 title.str("");
0185 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0186 hPullGP_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0187 title.str("");
0188 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0189 hPullGP_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0190
0191 title.str("");
0192 title << "PullGM_X_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0193 hPullGM_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0194 title.str("");
0195 title << "PullGM_Y_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0196 hPullGM_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0197 title.str("");
0198 title << "PullGM_Z_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0199 hPullGM_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0200
0201 title.str("");
0202 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_tr_stereo";
0203 hPullGP_X_tr_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0204 title.str("");
0205 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_tr_stereo";
0206 hPullGP_Y_tr_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0207 title.str("");
0208 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_tr_stereo";
0209 hPullGP_Z_tr_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0210
0211 title.str("");
0212 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_rs_stereo";
0213 hPullGP_X_rs_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0214 title.str("");
0215 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_rs_stereo";
0216 hPullGP_Y_rs_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0217 title.str("");
0218 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_rs_stereo";
0219 hPullGP_Z_rs_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0220 }
0221 }
0222 hTotChi2Increment = new TH1F("TotChi2Increment", "TotChi2Increment", 1000, 0, 100);
0223 hTotChi2GoodHit = new TH1F("TotChi2GoodHit", "TotChi2GoodHit", 1000, 0, 100);
0224 hTotChi2BadHit = new TH1F("TotChi2BadHit", "TotChi2BadHit", 1000, 0, 100);
0225 hTotChi2DeltaHit = new TH1F("TotChi2DeltaHit", "TotChi2DeltaHit", 1000, 0, 100);
0226 hTotChi2NSharedHit = new TH1F("TotChi2NSharedHit", "TotChi2NSharedHit", 1000, 0, 100);
0227 hTotChi2SharedHit = new TH1F("TotChi2SharedHit", "TotChi2SharedHit", 1000, 0, 100);
0228 hProcess_vs_Chi2 = new TH2F("Process_vs_Chi2", "Process_vs_Chi2", 1000, 0, 100, 17, -0.5, 16.5);
0229 hClsize_vs_Chi2 = new TH2F("Clsize_vs_Chi2", "Clsize_vs_Chi2", 1000, 0, 100, 17, -0.5, 16.5);
0230 hPixClsize_vs_Chi2 = new TH2F("PixClsize_vs_Chi2", "PixClsize_vs_Chi2", 1000, 0, 100, 17, -0.5, 16.5);
0231 hPrjClsize_vs_Chi2 = new TH2F("PrjClsize_vs_Chi2", "PrjClsize_vs_Chi2", 1000, 0, 100, 17, -0.5, 16.5);
0232 hSt1Clsize_vs_Chi2 = new TH2F("St1Clsize_vs_Chi2", "St1Clsize_vs_Chi2", 1000, 0, 100, 17, -0.5, 16.5);
0233 hSt2Clsize_vs_Chi2 = new TH2F("St2Clsize_vs_Chi2", "St2Clsize_vs_Chi2", 1000, 0, 100, 17, -0.5, 16.5);
0234 hGoodHit_vs_Chi2 = new TH2F("GoodHit_vs_Chi2", "GoodHit_vs_Chi2", 10000, 0, 1000, 2, -0.5, 1.5);
0235 hClusterSize = new TH1F("ClusterSize", "ClusterSize", 40, -0.5, 39.5);
0236 hPixClusterSize = new TH1F("PixClusterSize", "PixClusterSize", 40, -0.5, 39.5);
0237 hPrjClusterSize = new TH1F("PrjClusterSize", "PrjClusterSize", 40, -0.5, 39.5);
0238 hSt1ClusterSize = new TH1F("St1ClusterSize", "St1ClusterSize", 40, -0.5, 39.5);
0239 hSt2ClusterSize = new TH1F("St2ClusterSize", "St2ClusterSize", 40, -0.5, 39.5);
0240 hSimHitVecSize = new TH1F("hSimHitVecSize", "hSimHitVecSize", 40, -0.5, 39.5);
0241 hPixSimHitVecSize = new TH1F("PixSimHitVecSize", "PixSimHitVecSize", 40, -0.5, 39.5);
0242 hPrjSimHitVecSize = new TH1F("PrjSimHitVecSize", "PrjSimHitVecSize", 40, -0.5, 39.5);
0243 hSt1SimHitVecSize = new TH1F("St1SimHitVecSize", "St1SimHitVecSize", 40, -0.5, 39.5);
0244 hSt2SimHitVecSize = new TH1F("St2SimHitVecSize", "St2SimHitVecSize", 40, -0.5, 39.5);
0245 goodbadmerged = new TH1F("goodbadmerged", "goodbadmerged", 5, 0.5, 5.5);
0246 energyLossRatio = new TH1F("energyLossRatio", "energyLossRatio", 100, 0, 1);
0247 mergedPull = new TH1F("mergedPull", "mergedPull", 200, 0, 20);
0248 probXgood = new TH1F("probXgood", "probXgood", 110, 0, 1.1);
0249 probXbad = new TH1F("probXbad", "probXbad", 110, 0, 1.1);
0250 probXdelta = new TH1F("probXdelta", "probXdelta", 110, 0, 1.1);
0251 probXshared = new TH1F("probXshared", "probXshared", 110, 0, 1.1);
0252 probXnoshare = new TH1F("probXnoshare", "probXnoshare", 110, 0, 1.1);
0253 probYgood = new TH1F("probYgood", "probYgood", 110, 0, 1.1);
0254 probYbad = new TH1F("probYbad", "probYbad", 110, 0, 1.1);
0255 probYdelta = new TH1F("probYdelta", "probYdelta", 110, 0, 1.1);
0256 probYshared = new TH1F("probYshared", "probYshared", 110, 0, 1.1);
0257 probYnoshare = new TH1F("probYnoshare", "probYnoshare", 110, 0, 1.1);
0258 }
0259
0260 void TestTrackHits::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0261
0262 edm::ESHandle<TrackerTopology> tTopo = iSetup.getHandle(theTopoToken);
0263
0264 LogDebug("TestTrackHits") << "new event";
0265 iEvent.getByLabel(srcName, trajCollectionHandle);
0266 iEvent.getByLabel(srcName, trackCollectionHandle);
0267 iEvent.getByLabel(srcName, trajTrackAssociationCollectionHandle);
0268 iEvent.getByLabel(tpName, trackingParticleCollectionHandle);
0269 edm::Handle<reco::BeamSpot> beamSpot;
0270 iEvent.getByLabel("offlineBeamSpot", beamSpot);
0271
0272 LogTrace("TestTrackHits") << "Tr collection size=" << trackCollectionHandle->size();
0273 LogTrace("TestTrackHits") << "TP collection size=" << trackingParticleCollectionHandle->size();
0274
0275 iEvent.getByLabel("trackAssociatorByHits", trackAssociator);
0276
0277 TrackerHitAssociator hitAssociator(iEvent, trackerHitAssociatorConfig_);
0278
0279 reco::RecoToSimCollection recSimColl =
0280 trackAssociator->associateRecoToSim(trackCollectionHandle, trackingParticleCollectionHandle);
0281
0282 TrajectoryStateCombiner combiner;
0283
0284 int evtHits = 0;
0285 int i = 0;
0286 int yy = 0;
0287 int yyy = 0;
0288 for (std::vector<Trajectory>::const_iterator it = trajCollectionHandle->begin(); it != trajCollectionHandle->end();
0289 it++) {
0290 LogTrace("TestTrackHits") << "\n*****************new trajectory********************";
0291 double tchi2 = 0;
0292
0293 std::vector<TrajectoryMeasurement> tmColl = it->measurements();
0294
0295 edm::Ref<std::vector<Trajectory> > traj(trajCollectionHandle, i);
0296 reco::TrackRef tmptrack = (*trajTrackAssociationCollectionHandle.product())[traj];
0297 edm::RefToBase<reco::Track> track(tmptrack);
0298
0299
0300
0301
0302
0303
0304 std::vector<std::pair<TrackingParticleRef, double> > tP;
0305 if (recSimColl.find(track) != recSimColl.end()) {
0306 tP = recSimColl[track];
0307 if (!tP.empty()) {
0308 edm::LogVerbatim("TestTrackHits") << "reco::Track #" << ++yyy << " with pt=" << track->pt()
0309 << " associated with quality:" << tP.begin()->second << " good track #"
0310 << ++yy << " has hits:" << track->numberOfValidHits() << "\n";
0311 }
0312 } else {
0313 edm::LogVerbatim("TestTrackHits") << "reco::Track #" << ++yyy << " with pt=" << track->pt()
0314 << " NOT associated to any TrackingParticle"
0315 << "\n";
0316 i++;
0317 continue;
0318 }
0319
0320
0321
0322
0323
0324
0325
0326
0327
0328
0329 TrackingParticleRef tp = tP.begin()->first;
0330 LogTrace("TestTrackHits") << "a tp is associated with fraction=" << tP.begin()->second;
0331
0332 std::vector<unsigned int> tpids;
0333 for (TrackingParticle::g4t_iterator g4T = tp->g4Track_begin(); g4T != tp->g4Track_end(); ++g4T) {
0334 LogTrace("TestTrackHits") << "tp id=" << g4T->trackId();
0335 tpids.push_back(g4T->trackId());
0336 }
0337
0338
0339 int pp = 0;
0340 for (std::vector<TrajectoryMeasurement>::iterator tm = tmColl.begin(); tm != tmColl.end(); ++tm) {
0341 tchi2 += tm->estimate();
0342
0343 LogTrace("TestTrackHits") << "+++++++++++++++++new hit+++++++++++++++++";
0344 CTTRHp rhit = tm->recHit();
0345
0346
0347 TSOS state = combiner(tm->backwardPredictedState(), tm->forwardPredictedState());
0348
0349 if (rhit->isValid() == 0 && rhit->det() != nullptr)
0350 continue;
0351 evtHits++;
0352 LogTrace("TestTrackHits") << "valid hit #" << ++pp << "of hits=" << track->numberOfValidHits();
0353
0354 int subdetId = rhit->det()->geographicalId().subdetId();
0355 DetId id = rhit->det()->geographicalId();
0356 int layerId = tTopo->layer(id);
0357 LogTrace("TestTrackHits") << "subdetId=" << subdetId << " layerId=" << layerId;
0358
0359 const Surface* surf = rhit->surface();
0360 if (surf == nullptr)
0361 continue;
0362
0363 double energyLoss_ = 0.;
0364 unsigned int monoId = 0;
0365 std::vector<double> energyLossM;
0366 std::vector<double> energyLossS;
0367 std::vector<PSimHit> assSimHits = hitAssociator.associateHit(*(rhit)->hit());
0368 unsigned int simhitvecsize = assSimHits.size();
0369 if (simhitvecsize == 0)
0370 continue;
0371 PSimHit shit;
0372 std::vector<unsigned int> trackIds;
0373 energyLossS.clear();
0374 energyLossM.clear();
0375 for (std::vector<PSimHit>::const_iterator m = assSimHits.begin(); m < assSimHits.end(); m++) {
0376 unsigned int tId = m->trackId();
0377 if (find(trackIds.begin(), trackIds.end(), tId) == trackIds.end())
0378 trackIds.push_back(tId);
0379 if (m->energyLoss() > energyLoss_) {
0380 shit = *m;
0381 energyLoss_ = m->energyLoss();
0382 }
0383 if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
0384 if (monoId == 0)
0385 monoId = m->detUnitId();
0386 if (monoId == m->detUnitId()) {
0387 energyLossM.push_back(m->energyLoss());
0388 } else {
0389 energyLossS.push_back(m->energyLoss());
0390 }
0391
0392 } else {
0393 energyLossM.push_back(m->energyLoss());
0394 }
0395 }
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406
0407
0408
0409
0410 double chi2increment = tm->estimate();
0411 LogTrace("TestTrackHits") << "tm->estimate()=" << tm->estimate();
0412 title.str("");
0413 title << "Chi2Increment_" << subdetId << "-" << layerId;
0414 hChi2Increment[title.str()]->Fill(chi2increment);
0415 title.str("");
0416 title << "Chi2IncrementVsEta_" << subdetId << "-" << layerId;
0417 hChi2IncrementVsEta[title.str()]->Fill(track->eta(), chi2increment);
0418 hTotChi2Increment->Fill(chi2increment);
0419 hProcess_vs_Chi2->Fill(chi2increment, shit.processType());
0420
0421 int clustersize = 0;
0422 bool mergedhit = false;
0423 if (dynamic_cast<const SiPixelRecHit*>(rhit->hit())) {
0424 clustersize = ((const SiPixelRecHit*)(rhit->hit()))->cluster()->size();
0425 hPixClsize_vs_Chi2->Fill(chi2increment, clustersize);
0426 hPixClusterSize->Fill(clustersize);
0427 hPixSimHitVecSize->Fill(simhitvecsize);
0428 if (simhitvecsize > 1)
0429 mergedhit = true;
0430 } else if (dynamic_cast<const SiStripRecHit2D*>(rhit->hit())) {
0431 clustersize = ((const SiStripRecHit2D*)(rhit->hit()))->cluster()->amplitudes().size();
0432 hSt1Clsize_vs_Chi2->Fill(chi2increment, clustersize);
0433 hSt1ClusterSize->Fill(clustersize);
0434 hSt1SimHitVecSize->Fill(simhitvecsize);
0435 if (simhitvecsize > 1)
0436 mergedhit = true;
0437 } else if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
0438 int clsize1 = ((const SiStripMatchedRecHit2D*)(rhit->hit()))->monoCluster().amplitudes().size();
0439 int clsize2 = ((const SiStripMatchedRecHit2D*)(rhit->hit()))->stereoCluster().amplitudes().size();
0440 if (clsize1 > clsize2)
0441 clustersize = clsize1;
0442 else
0443 clustersize = clsize2;
0444 hSt2Clsize_vs_Chi2->Fill(chi2increment, clustersize);
0445 hSt2ClusterSize->Fill(clustersize);
0446 hSt2SimHitVecSize->Fill(simhitvecsize);
0447 if (simhitvecsize > 2)
0448 mergedhit = true;
0449 } else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(rhit->hit())) {
0450 clustersize = ((const ProjectedSiStripRecHit2D*)(rhit->hit()))->originalHit().cluster()->amplitudes().size();
0451 hPrjClsize_vs_Chi2->Fill(chi2increment, clustersize);
0452 hPrjClusterSize->Fill(clustersize);
0453 hPrjSimHitVecSize->Fill(simhitvecsize);
0454 if (simhitvecsize > 1)
0455 mergedhit = true;
0456 }
0457 hClsize_vs_Chi2->Fill(chi2increment, clustersize);
0458 hClusterSize->Fill(clustersize);
0459 hSimHitVecSize->Fill(simhitvecsize);
0460
0461
0462
0463
0464
0465
0466 std::vector<SimHitIdpr> simTrackIds = hitAssociator.associateHitId(*(rhit)->hit());
0467 bool goodhit = false;
0468 for (size_t j = 0; j < simTrackIds.size(); j++) {
0469 LogTrace("TestTrackHits") << "hit id=" << simTrackIds[j].first;
0470 for (size_t jj = 0; jj < tpids.size(); jj++) {
0471 if (simTrackIds[j].first == tpids[jj])
0472 goodhit = true;
0473 break;
0474 }
0475 }
0476 bool shared = true;
0477 bool ioniOnly = true;
0478 const SiPixelRecHit* pix = dynamic_cast<const SiPixelRecHit*>(rhit->hit());
0479 if (goodhit) {
0480 if (energyLossM.size() > 1 && energyLossS.size() <= 1) {
0481 sort(energyLossM.begin(), energyLossM.end(), greater<double>());
0482 energyLossRatio->Fill(energyLossM[1] / energyLossM[0]);
0483 } else if (energyLossS.size() > 1 && energyLossM.size() <= 1) {
0484 sort(energyLossS.begin(), energyLossS.end(), greater<double>());
0485 energyLossRatio->Fill(energyLossS[1] / energyLossS[0]);
0486 } else if (energyLossS.size() > 1 && energyLossM.size() > 1) {
0487 sort(energyLossM.begin(), energyLossM.end(), greater<double>());
0488 sort(energyLossS.begin(), energyLossS.end(), greater<double>());
0489 energyLossM[1] / energyLossM[0] > energyLossS[1] / energyLossS[0]
0490 ? energyLossRatio->Fill(energyLossM[1] / energyLossM[0])
0491 : energyLossRatio->Fill(energyLossS[1] / energyLossS[0]);
0492 }
0493 if (mergedhit) {
0494
0495 LogVerbatim("TestTrackHits") << "MERGED HIT" << std::endl;
0496 unsigned int idc = 0;
0497 for (size_t jj = 0; jj < tpids.size(); jj++) {
0498 idc += std::count(trackIds.begin(), trackIds.end(), tpids[jj]);
0499 }
0500 if (idc == trackIds.size()) {
0501 shared = false;
0502 }
0503 for (std::vector<PSimHit>::const_iterator m = assSimHits.begin() + 1; m < assSimHits.end(); m++) {
0504 if ((m->processType() != 7 && m->processType() != 8 && m->processType() != 9) &&
0505 abs(m->particleType()) != 11) {
0506 ioniOnly = false;
0507 break;
0508 }
0509 }
0510 if (ioniOnly && !shared) {
0511 title.str("");
0512 title << "Chi2DeltaHit_" << subdetId << "-" << layerId;
0513 hChi2DeltaHit[title.str()]->Fill(chi2increment);
0514 hTotChi2DeltaHit->Fill(chi2increment);
0515 if (pix) {
0516 probXdelta->Fill(pix->probabilityX());
0517 probYdelta->Fill(pix->probabilityY());
0518 }
0519 } else if (!ioniOnly && !shared) {
0520 title.str("");
0521 title << "Chi2NSharedHit_" << subdetId << "-" << layerId;
0522 hChi2NSharedHit[title.str()]->Fill(chi2increment);
0523 hTotChi2NSharedHit->Fill(chi2increment);
0524 if (pix) {
0525 probXnoshare->Fill(pix->probabilityX());
0526 probYnoshare->Fill(pix->probabilityY());
0527 }
0528 } else {
0529 title.str("");
0530 title << "Chi2SharedHit_" << subdetId << "-" << layerId;
0531 hChi2SharedHit[title.str()]->Fill(chi2increment);
0532 hTotChi2SharedHit->Fill(chi2increment);
0533 if (pix) {
0534 probXshared->Fill(pix->probabilityX());
0535 probYshared->Fill(pix->probabilityY());
0536 }
0537 }
0538
0539 for (std::vector<PSimHit>::const_iterator m = assSimHits.begin(); m < assSimHits.end(); m++) {
0540 unsigned int tId = m->trackId();
0541 LogVerbatim("TestTrackHits") << "component with id=" << tId << " eLoss=" << m->energyLoss()
0542 << " pType=" << m->processType();
0543 if (find(tpids.begin(), tpids.end(), tId) == tpids.end())
0544 continue;
0545 if (m->processType() == 2) {
0546 GlobalPoint gpr = rhit->globalPosition();
0547 AlgebraicSymMatrix33 ger = rhit->globalPositionError().matrix();
0548 GlobalPoint gps = surf->toGlobal(m->localPosition());
0549 LogVerbatim("TestTrackHits") << gpr << " " << gps << " " << ger;
0550 ROOT::Math::SVector<double, 3> delta;
0551 delta[0] = gpr.x() - gps.x();
0552 delta[1] = gpr.y() - gps.y();
0553 delta[2] = gpr.z() - gps.z();
0554 LogVerbatim("TestTrackHits") << delta << " " << ger;
0555 double mpull = sqrt(delta[0] * delta[0] / ger[0][0] + delta[1] * delta[1] / ger[1][1] +
0556 delta[2] * delta[2] / ger[2][2]);
0557 LogVerbatim("TestTrackHits") << "hit pull=" << mpull;
0558 mergedPull->Fill(mpull);
0559 break;
0560 }
0561 }
0562 } else {
0563 LogVerbatim("TestTrackHits") << "good hit";
0564 title.str("");
0565 title << "Chi2GoodHit_" << subdetId << "-" << layerId;
0566 hChi2GoodHit[title.str()]->Fill(chi2increment);
0567 hTotChi2GoodHit->Fill(chi2increment);
0568 if (pix) {
0569 probXgood->Fill(pix->probabilityX());
0570 probYgood->Fill(pix->probabilityY());
0571 }
0572 }
0573 } else {
0574 LogVerbatim("TestTrackHits") << "bad hit";
0575 title.str("");
0576 title << "Chi2BadHit_" << subdetId << "-" << layerId;
0577 hChi2BadHit[title.str()]->Fill(chi2increment);
0578 hTotChi2BadHit->Fill(chi2increment);
0579 goodbadmerged->Fill(2);
0580 if (pix) {
0581 probXbad->Fill(pix->probabilityX());
0582 probYbad->Fill(pix->probabilityY());
0583 }
0584 }
0585 hGoodHit_vs_Chi2->Fill(chi2increment, goodhit);
0586
0587 LocalVector shitLMom;
0588 LocalPoint shitLPos;
0589 if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
0590 if (simhitvecsize > 2 && goodhit) {
0591 if (ioniOnly && !shared)
0592 goodbadmerged->Fill(3);
0593 else if (!ioniOnly && !shared)
0594 goodbadmerged->Fill(4);
0595 else
0596 goodbadmerged->Fill(5);
0597 } else if (goodhit)
0598 goodbadmerged->Fill(1);
0599 double rechitmatchedx = rhit->localPosition().x();
0600 double rechitmatchedy = rhit->localPosition().y();
0601 double mindist = 999999;
0602 float distx, disty;
0603 std::pair<LocalPoint, LocalVector> closestPair;
0604 const StripGeomDetUnit* stripDet = (StripGeomDetUnit*)((const GluedGeomDet*)(rhit)->det())->stereoDet();
0605 const BoundPlane& plane = (rhit)->det()->surface();
0606 for (std::vector<PSimHit>::const_iterator m = assSimHits.begin(); m < assSimHits.end(); m++) {
0607
0608 std::pair<LocalPoint, LocalVector> hitPair = projectHit((*m), stripDet, plane);
0609 distx = fabs(rechitmatchedx - hitPair.first.x());
0610 disty = fabs(rechitmatchedy - hitPair.first.y());
0611 double dist = distx * distx + disty * disty;
0612 if (sqrt(dist) < mindist) {
0613 mindist = dist;
0614 closestPair = hitPair;
0615 }
0616 }
0617 shitLPos = closestPair.first;
0618 shitLMom = closestPair.second;
0619 } else {
0620 if (simhitvecsize > 1 && goodhit) {
0621 if (ioniOnly && !shared)
0622 goodbadmerged->Fill(3);
0623 else if (!ioniOnly && !shared)
0624 goodbadmerged->Fill(4);
0625 else
0626 goodbadmerged->Fill(5);
0627 } else if (goodhit)
0628 goodbadmerged->Fill(1);
0629 shitLPos = shit.localPosition();
0630 shitLMom = shit.momentumAtEntry();
0631 }
0632 GlobalVector shitGMom = surf->toGlobal(shitLMom);
0633 GlobalPoint shitGPos = surf->toGlobal(shitLPos);
0634
0635 GlobalVector tsosGMom = state.globalMomentum();
0636 GlobalError tsosGMEr(state.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3, 3));
0637 GlobalPoint tsosGPos = state.globalPosition();
0638 GlobalError tsosGPEr = state.cartesianError().position();
0639
0640 GlobalPoint rhitGPos = (rhit)->globalPosition();
0641 GlobalError rhitGPEr = (rhit)->globalPositionError();
0642
0643 LogVerbatim("TestTrackHits") << "assSimHits.size()=" << assSimHits.size();
0644 LogVerbatim("TestTrackHits") << "tsos globalPos =" << tsosGPos;
0645 LogVerbatim("TestTrackHits") << "sim hit globalPos=" << shitGPos;
0646 LogVerbatim("TestTrackHits") << "rec hit globalPos=" << rhitGPos;
0647 LogVerbatim("TestTrackHits") << "geographicalId =" << rhit->det()->geographicalId().rawId();
0648 LogVerbatim("TestTrackHits") << "surface position =" << surf->position();
0649
0650 #if 0
0651 if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->geographicalId()="
0652 << rhit->detUnit()->geographicalId().rawId() ;
0653 LogTrace("TestTrackHits") << "rhit->det()->surface().position()="
0654 << rhit->det()->surface().position() ;
0655 if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().position()="
0656 << rhit->detUnit()->surface().position() ;
0657 LogTrace("TestTrackHits") << "rhit->det()->position()=" << rhit->det()->position() ;
0658 if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->position()="
0659 << rhit->detUnit()->position() ;
0660 LogTrace("TestTrackHits") << "rhit->det()->surface().bounds().length()="
0661 << rhit->det()->surface().bounds().length() ;
0662 if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().bounds().length()="
0663 << rhit->detUnit()->surface().bounds().length() ;
0664 LogTrace("TestTrackHits") << "rhit->det()->surface().bounds().width()="
0665 << rhit->det()->surface().bounds().width() ;
0666 if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().bounds().width()="
0667 << rhit->detUnit()->surface().bounds().width() ;
0668 LogTrace("TestTrackHits") << "rhit->det()->surface().bounds().thickness()="
0669 << rhit->det()->surface().bounds().thickness() ;
0670 if (rhit->detUnit()) LogTrace("TestTrackHits") << "rhit->detUnit()->surface().bounds().thickness()="
0671 << rhit->detUnit()->surface().bounds().thickness() ;
0672 #endif
0673
0674 double pullGPX_rs = (rhitGPos.x() - shitGPos.x()) / sqrt(rhitGPEr.cxx());
0675 double pullGPY_rs = (rhitGPos.y() - shitGPos.y()) / sqrt(rhitGPEr.cyy());
0676 double pullGPZ_rs = (rhitGPos.z() - shitGPos.z()) / sqrt(rhitGPEr.czz());
0677
0678
0679
0680
0681 LogTrace("TestTrackHits") << "rs";
0682
0683 title.str("");
0684 title << "PullGP_X_" << subdetId << "-" << layerId << "_rs";
0685 hPullGP_X_rs[title.str()]->Fill(pullGPX_rs);
0686 title.str("");
0687 title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs";
0688 hPullGP_Y_rs[title.str()]->Fill(pullGPY_rs);
0689 title.str("");
0690 title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs";
0691 hPullGP_Z_rs[title.str()]->Fill(pullGPZ_rs);
0692
0693 double pullGPX_tr = (tsosGPos.x() - rhitGPos.x()) / sqrt(tsosGPEr.cxx() + rhitGPEr.cxx());
0694 double pullGPY_tr = (tsosGPos.y() - rhitGPos.y()) / sqrt(tsosGPEr.cyy() + rhitGPEr.cyy());
0695 double pullGPZ_tr = (tsosGPos.z() - rhitGPos.z()) / sqrt(tsosGPEr.czz() + rhitGPEr.czz());
0696
0697
0698
0699
0700 LogTrace("TestTrackHits") << "tr";
0701
0702 title.str("");
0703 title << "PullGP_X_" << subdetId << "-" << layerId << "_tr";
0704 hPullGP_X_tr[title.str()]->Fill(pullGPX_tr);
0705 title.str("");
0706 title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr";
0707 hPullGP_Y_tr[title.str()]->Fill(pullGPY_tr);
0708 title.str("");
0709 title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr";
0710 hPullGP_Z_tr[title.str()]->Fill(pullGPZ_tr);
0711
0712 double pullGPX_ts = (tsosGPos.x() - shitGPos.x()) / sqrt(tsosGPEr.cxx());
0713 double pullGPY_ts = (tsosGPos.y() - shitGPos.y()) / sqrt(tsosGPEr.cyy());
0714 double pullGPZ_ts = (tsosGPos.z() - shitGPos.z()) / sqrt(tsosGPEr.czz());
0715
0716
0717
0718
0719 LogTrace("TestTrackHits") << "ts1";
0720
0721 title.str("");
0722 title << "PullGP_X_" << subdetId << "-" << layerId << "_ts";
0723 hPullGP_X_ts[title.str()]->Fill(pullGPX_ts);
0724 title.str("");
0725 title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts";
0726 hPullGP_Y_ts[title.str()]->Fill(pullGPY_ts);
0727 title.str("");
0728 title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts";
0729 hPullGP_Z_ts[title.str()]->Fill(pullGPZ_ts);
0730
0731 double pullGMX_ts = (tsosGMom.x() - shitGMom.x()) / sqrt(tsosGMEr.cxx());
0732 double pullGMY_ts = (tsosGMom.y() - shitGMom.y()) / sqrt(tsosGMEr.cyy());
0733 double pullGMZ_ts = (tsosGMom.z() - shitGMom.z()) / sqrt(tsosGMEr.czz());
0734
0735
0736
0737
0738 LogTrace("TestTrackHits") << "ts2";
0739
0740 title.str("");
0741 title << "PullGM_X_" << subdetId << "-" << layerId << "_ts";
0742 hPullGM_X_ts[title.str()]->Fill(pullGMX_ts);
0743 title.str("");
0744 title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts";
0745 hPullGM_Y_ts[title.str()]->Fill(pullGMY_ts);
0746 title.str("");
0747 title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts";
0748 hPullGM_Z_ts[title.str()]->Fill(pullGMZ_ts);
0749 if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
0750 Propagator* thePropagatorAnyDir = new PropagatorWithMaterial(anyDirection, 0.105, theMF.product(), 1.6);
0751
0752 LogTrace("TestTrackHits") << "MONO HIT";
0753 auto m = dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())->monoHit();
0754 CTTRHp tMonoHit = theBuilder->build(&m);
0755 if (tMonoHit == nullptr)
0756 continue;
0757 vector<PSimHit> assMonoSimHits = hitAssociator.associateHit(*tMonoHit->hit());
0758 if (assMonoSimHits.empty())
0759 continue;
0760 const PSimHit sMonoHit = *(assMonoSimHits.begin());
0761 const Surface* monoSurf = &(tMonoHit->det()->surface());
0762 if (monoSurf == nullptr)
0763 continue;
0764 TSOS monoState = thePropagatorAnyDir->propagate(state, *monoSurf);
0765 if (monoState.isValid() == 0)
0766 continue;
0767
0768 LocalVector monoShitLMom = sMonoHit.momentumAtEntry();
0769 GlobalVector monoShitGMom = monoSurf->toGlobal(monoShitLMom);
0770 LocalPoint monoShitLPos = sMonoHit.localPosition();
0771 GlobalPoint monoShitGPos = monoSurf->toGlobal(monoShitLPos);
0772
0773
0774
0775
0776 GlobalVector monoTsosGMom = monoState.globalMomentum();
0777 GlobalError monoTsosGMEr(monoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3, 3));
0778 GlobalPoint monoTsosGPos = monoState.globalPosition();
0779 GlobalError monoTsosGPEr = monoState.cartesianError().position();
0780
0781 GlobalPoint monoRhitGPos = tMonoHit->globalPosition();
0782 GlobalError monoRhitGPEr = tMonoHit->globalPositionError();
0783
0784 double pullGPX_rs_mono = (monoRhitGPos.x() - monoShitGPos.x()) / sqrt(monoRhitGPEr.cxx());
0785 double pullGPY_rs_mono = (monoRhitGPos.y() - monoShitGPos.y()) / sqrt(monoRhitGPEr.cyy());
0786 double pullGPZ_rs_mono = (monoRhitGPos.z() - monoShitGPos.z()) / sqrt(monoRhitGPEr.czz());
0787
0788
0789
0790
0791 MeasurementExtractor meMo(monoState);
0792 double chi2mono = computeChi2Increment(meMo, tMonoHit);
0793
0794 title.str("");
0795 title << "Chi2Increment_mono_" << subdetId << "-" << layerId;
0796 hChi2Increment_mono[title.str()]->Fill(chi2mono);
0797
0798 title.str("");
0799 title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_mono";
0800 hPullGP_X_rs_mono[title.str()]->Fill(pullGPX_rs_mono);
0801 title.str("");
0802 title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_mono";
0803 hPullGP_Y_rs_mono[title.str()]->Fill(pullGPY_rs_mono);
0804 title.str("");
0805 title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_mono";
0806 hPullGP_Z_rs_mono[title.str()]->Fill(pullGPZ_rs_mono);
0807
0808 double pullGPX_tr_mono = (monoTsosGPos.x() - monoRhitGPos.x()) / sqrt(monoTsosGPEr.cxx() + monoRhitGPEr.cxx());
0809 double pullGPY_tr_mono = (monoTsosGPos.y() - monoRhitGPos.y()) / sqrt(monoTsosGPEr.cyy() + monoRhitGPEr.cyy());
0810 double pullGPZ_tr_mono = (monoTsosGPos.z() - monoRhitGPos.z()) / sqrt(monoTsosGPEr.czz() + monoRhitGPEr.czz());
0811
0812
0813
0814
0815 title.str("");
0816 title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_mono";
0817 hPullGP_X_tr_mono[title.str()]->Fill(pullGPX_tr_mono);
0818 title.str("");
0819 title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_mono";
0820 hPullGP_Y_tr_mono[title.str()]->Fill(pullGPY_tr_mono);
0821 title.str("");
0822 title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_mono";
0823 hPullGP_Z_tr_mono[title.str()]->Fill(pullGPZ_tr_mono);
0824
0825 double pullGPX_ts_mono = (monoTsosGPos.x() - monoShitGPos.x()) / sqrt(monoTsosGPEr.cxx());
0826 double pullGPY_ts_mono = (monoTsosGPos.y() - monoShitGPos.y()) / sqrt(monoTsosGPEr.cyy());
0827 double pullGPZ_ts_mono = (monoTsosGPos.z() - monoShitGPos.z()) / sqrt(monoTsosGPEr.czz());
0828
0829
0830
0831
0832 title.str("");
0833 title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_mono";
0834 hPullGP_X_ts_mono[title.str()]->Fill(pullGPX_ts_mono);
0835 title.str("");
0836 title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_mono";
0837 hPullGP_Y_ts_mono[title.str()]->Fill(pullGPY_ts_mono);
0838 title.str("");
0839 title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_mono";
0840 hPullGP_Z_ts_mono[title.str()]->Fill(pullGPZ_ts_mono);
0841
0842 double pullGMX_ts_mono = (monoTsosGMom.x() - monoShitGMom.x()) / sqrt(monoTsosGMEr.cxx());
0843 double pullGMY_ts_mono = (monoTsosGMom.y() - monoShitGMom.y()) / sqrt(monoTsosGMEr.cyy());
0844 double pullGMZ_ts_mono = (monoTsosGMom.z() - monoShitGMom.z()) / sqrt(monoTsosGMEr.czz());
0845
0846
0847
0848
0849 title.str("");
0850 title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_mono";
0851 hPullGM_X_ts_mono[title.str()]->Fill(pullGMX_ts_mono);
0852 title.str("");
0853 title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_mono";
0854 hPullGM_Y_ts_mono[title.str()]->Fill(pullGMY_ts_mono);
0855 title.str("");
0856 title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_mono";
0857 hPullGM_Z_ts_mono[title.str()]->Fill(pullGMZ_ts_mono);
0858
0859
0860 LogTrace("TestTrackHits") << "STEREO HIT";
0861 auto s = dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())->stereoHit();
0862 CTTRHp tStereoHit = theBuilder->build(&s);
0863 if (tStereoHit == nullptr)
0864 continue;
0865 vector<PSimHit> assStereoSimHits = hitAssociator.associateHit(*tStereoHit->hit());
0866 if (assStereoSimHits.empty())
0867 continue;
0868 const PSimHit sStereoHit = *(assStereoSimHits.begin());
0869 const Surface* stereoSurf = &(tStereoHit->det()->surface());
0870 if (stereoSurf == nullptr)
0871 continue;
0872 TSOS stereoState = thePropagatorAnyDir->propagate(state, *stereoSurf);
0873 if (stereoState.isValid() == 0)
0874 continue;
0875
0876 LocalVector stereoShitLMom = sStereoHit.momentumAtEntry();
0877 GlobalVector stereoShitGMom = stereoSurf->toGlobal(stereoShitLMom);
0878 LocalPoint stereoShitLPos = sStereoHit.localPosition();
0879 GlobalPoint stereoShitGPos = stereoSurf->toGlobal(stereoShitLPos);
0880
0881
0882
0883
0884 GlobalVector stereoTsosGMom = stereoState.globalMomentum();
0885 GlobalError stereoTsosGMEr(stereoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3, 3));
0886 GlobalPoint stereoTsosGPos = stereoState.globalPosition();
0887 GlobalError stereoTsosGPEr = stereoState.cartesianError().position();
0888
0889 GlobalPoint stereoRhitGPos = tStereoHit->globalPosition();
0890 GlobalError stereoRhitGPEr = tStereoHit->globalPositionError();
0891
0892 MeasurementExtractor meSt(stereoState);
0893 double chi2stereo = computeChi2Increment(meSt, tStereoHit);
0894
0895 title.str("");
0896 title << "Chi2Increment_stereo_" << subdetId << "-" << layerId;
0897 hChi2Increment_stereo[title.str()]->Fill(chi2stereo);
0898
0899 double pullGPX_rs_stereo = (stereoRhitGPos.x() - stereoShitGPos.x()) / sqrt(stereoRhitGPEr.cxx());
0900 double pullGPY_rs_stereo = (stereoRhitGPos.y() - stereoShitGPos.y()) / sqrt(stereoRhitGPEr.cyy());
0901 double pullGPZ_rs_stereo = (stereoRhitGPos.z() - stereoShitGPos.z()) / sqrt(stereoRhitGPEr.czz());
0902
0903
0904
0905
0906 title.str("");
0907 title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_stereo";
0908 hPullGP_X_rs_stereo[title.str()]->Fill(pullGPX_rs_stereo);
0909 title.str("");
0910 title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_stereo";
0911 hPullGP_Y_rs_stereo[title.str()]->Fill(pullGPY_rs_stereo);
0912 title.str("");
0913 title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_stereo";
0914 hPullGP_Z_rs_stereo[title.str()]->Fill(pullGPZ_rs_stereo);
0915
0916 double pullGPX_tr_stereo =
0917 (stereoTsosGPos.x() - stereoRhitGPos.x()) / sqrt(stereoTsosGPEr.cxx() + stereoRhitGPEr.cxx());
0918 double pullGPY_tr_stereo =
0919 (stereoTsosGPos.y() - stereoRhitGPos.y()) / sqrt(stereoTsosGPEr.cyy() + stereoRhitGPEr.cyy());
0920 double pullGPZ_tr_stereo =
0921 (stereoTsosGPos.z() - stereoRhitGPos.z()) / sqrt(stereoTsosGPEr.czz() + stereoRhitGPEr.czz());
0922
0923
0924
0925
0926 title.str("");
0927 title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_stereo";
0928 hPullGP_X_tr_stereo[title.str()]->Fill(pullGPX_tr_stereo);
0929 title.str("");
0930 title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_stereo";
0931 hPullGP_Y_tr_stereo[title.str()]->Fill(pullGPY_tr_stereo);
0932 title.str("");
0933 title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_stereo";
0934 hPullGP_Z_tr_stereo[title.str()]->Fill(pullGPZ_tr_stereo);
0935
0936 double pullGPX_ts_stereo = (stereoTsosGPos.x() - stereoShitGPos.x()) / sqrt(stereoTsosGPEr.cxx());
0937 double pullGPY_ts_stereo = (stereoTsosGPos.y() - stereoShitGPos.y()) / sqrt(stereoTsosGPEr.cyy());
0938 double pullGPZ_ts_stereo = (stereoTsosGPos.z() - stereoShitGPos.z()) / sqrt(stereoTsosGPEr.czz());
0939
0940
0941
0942
0943 title.str("");
0944 title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_stereo";
0945 hPullGP_X_ts_stereo[title.str()]->Fill(pullGPX_ts_stereo);
0946 title.str("");
0947 title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_stereo";
0948 hPullGP_Y_ts_stereo[title.str()]->Fill(pullGPY_ts_stereo);
0949 title.str("");
0950 title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_stereo";
0951 hPullGP_Z_ts_stereo[title.str()]->Fill(pullGPZ_ts_stereo);
0952
0953 double pullGMX_ts_stereo = (stereoTsosGMom.x() - stereoShitGMom.x()) / sqrt(stereoTsosGMEr.cxx());
0954 double pullGMY_ts_stereo = (stereoTsosGMom.y() - stereoShitGMom.y()) / sqrt(stereoTsosGMEr.cyy());
0955 double pullGMZ_ts_stereo = (stereoTsosGMom.z() - stereoShitGMom.z()) / sqrt(stereoTsosGMEr.czz());
0956
0957
0958
0959
0960 title.str("");
0961 title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_stereo";
0962 hPullGM_X_ts_stereo[title.str()]->Fill(pullGMX_ts_stereo);
0963 title.str("");
0964 title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_stereo";
0965 hPullGM_Y_ts_stereo[title.str()]->Fill(pullGMY_ts_stereo);
0966 title.str("");
0967 title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_stereo";
0968 hPullGM_Z_ts_stereo[title.str()]->Fill(pullGMZ_ts_stereo);
0969 }
0970 }
0971 LogTrace("TestTrackHits") << "traj chi2=" << tchi2;
0972 LogTrace("TestTrackHits") << "track chi2=" << track->chi2();
0973 i++;
0974 }
0975 LogTrace("TestTrackHits") << "end of event: processd hits=" << evtHits;
0976 }
0977
0978 void TestTrackHits::endJob() {
0979
0980 TDirectory* chi2d = file->mkdir("Chi2Increment");
0981
0982 TDirectory* gp_ts = file->mkdir("GP_TSOS-SimHit");
0983 TDirectory* gm_ts = file->mkdir("GM_TSOS-SimHit");
0984 TDirectory* gp_tr = file->mkdir("GP_TSOS-RecHit");
0985 TDirectory* gp_rs = file->mkdir("GP_RecHit-SimHit");
0986
0987 TDirectory* gp_tsx = gp_ts->mkdir("X");
0988 TDirectory* gp_tsy = gp_ts->mkdir("Y");
0989 TDirectory* gp_tsz = gp_ts->mkdir("Z");
0990 TDirectory* gm_tsx = gm_ts->mkdir("X");
0991 TDirectory* gm_tsy = gm_ts->mkdir("Y");
0992 TDirectory* gm_tsz = gm_ts->mkdir("Z");
0993 TDirectory* gp_trx = gp_tr->mkdir("X");
0994 TDirectory* gp_try = gp_tr->mkdir("Y");
0995 TDirectory* gp_trz = gp_tr->mkdir("Z");
0996 TDirectory* gp_rsx = gp_rs->mkdir("X");
0997 TDirectory* gp_rsy = gp_rs->mkdir("Y");
0998 TDirectory* gp_rsz = gp_rs->mkdir("Z");
0999
1000 TDirectory* gp_tsx_mono = gp_ts->mkdir("MONOX");
1001 TDirectory* gp_tsy_mono = gp_ts->mkdir("MONOY");
1002 TDirectory* gp_tsz_mono = gp_ts->mkdir("MONOZ");
1003 TDirectory* gm_tsx_mono = gm_ts->mkdir("MONOX");
1004 TDirectory* gm_tsy_mono = gm_ts->mkdir("MONOY");
1005 TDirectory* gm_tsz_mono = gm_ts->mkdir("MONOZ");
1006 TDirectory* gp_trx_mono = gp_tr->mkdir("MONOX");
1007 TDirectory* gp_try_mono = gp_tr->mkdir("MONOY");
1008 TDirectory* gp_trz_mono = gp_tr->mkdir("MONOZ");
1009 TDirectory* gp_rsx_mono = gp_rs->mkdir("MONOX");
1010 TDirectory* gp_rsy_mono = gp_rs->mkdir("MONOY");
1011 TDirectory* gp_rsz_mono = gp_rs->mkdir("MONOZ");
1012
1013 TDirectory* gp_tsx_stereo = gp_ts->mkdir("STEREOX");
1014 TDirectory* gp_tsy_stereo = gp_ts->mkdir("STEREOY");
1015 TDirectory* gp_tsz_stereo = gp_ts->mkdir("STEREOZ");
1016 TDirectory* gm_tsx_stereo = gm_ts->mkdir("STEREOX");
1017 TDirectory* gm_tsy_stereo = gm_ts->mkdir("STEREOY");
1018 TDirectory* gm_tsz_stereo = gm_ts->mkdir("STEREOZ");
1019 TDirectory* gp_trx_stereo = gp_tr->mkdir("STEREOX");
1020 TDirectory* gp_try_stereo = gp_tr->mkdir("STEREOY");
1021 TDirectory* gp_trz_stereo = gp_tr->mkdir("STEREOZ");
1022 TDirectory* gp_rsx_stereo = gp_rs->mkdir("STEREOX");
1023 TDirectory* gp_rsy_stereo = gp_rs->mkdir("STEREOY");
1024 TDirectory* gp_rsz_stereo = gp_rs->mkdir("STEREOZ");
1025
1026 chi2d->cd();
1027 hTotChi2Increment->Write();
1028 hTotChi2GoodHit->Write();
1029 hTotChi2BadHit->Write();
1030 hTotChi2DeltaHit->Write();
1031 hTotChi2NSharedHit->Write();
1032 hTotChi2SharedHit->Write();
1033 hProcess_vs_Chi2->Write();
1034 hClsize_vs_Chi2->Write();
1035 hPixClsize_vs_Chi2->Write();
1036 hPrjClsize_vs_Chi2->Write();
1037 hSt1Clsize_vs_Chi2->Write();
1038 hSt2Clsize_vs_Chi2->Write();
1039 hGoodHit_vs_Chi2->Write();
1040 hClusterSize->Write();
1041 hPixClusterSize->Write();
1042 hPrjClusterSize->Write();
1043 hSt1ClusterSize->Write();
1044 hSt2ClusterSize->Write();
1045 hSimHitVecSize->Write();
1046 hPixSimHitVecSize->Write();
1047 hPrjSimHitVecSize->Write();
1048 hSt1SimHitVecSize->Write();
1049 hSt2SimHitVecSize->Write();
1050 goodbadmerged->Write();
1051 energyLossRatio->Write();
1052 mergedPull->Write();
1053 probXgood->Write();
1054 probXbad->Write();
1055 probXdelta->Write();
1056 probXshared->Write();
1057 probXnoshare->Write();
1058 probYgood->Write();
1059 probYbad->Write();
1060 probYdelta->Write();
1061 probYshared->Write();
1062 probYnoshare->Write();
1063 for (int i = 0; i != 6; i++)
1064 for (int j = 0; j != 9; j++) {
1065 if (i == 0 && j > 2)
1066 break;
1067 if (i == 1 && j > 1)
1068 break;
1069 if (i == 2 && j > 3)
1070 break;
1071 if (i == 3 && j > 2)
1072 break;
1073 if (i == 4 && j > 5)
1074 break;
1075 if (i == 5 && j > 8)
1076 break;
1077 chi2d->cd();
1078 title.str("");
1079 title << "Chi2Increment_" << i + 1 << "-" << j + 1;
1080 hChi2Increment[title.str()]->Write();
1081 title.str("");
1082 title << "Chi2IncrementVsEta_" << i + 1 << "-" << j + 1;
1083 hChi2IncrementVsEta[title.str()]->Write();
1084 title.str("");
1085 title << "Chi2GoodHit_" << i + 1 << "-" << j + 1;
1086 hChi2GoodHit[title.str()]->Write();
1087 title.str("");
1088 title << "Chi2BadHit_" << i + 1 << "-" << j + 1;
1089 hChi2BadHit[title.str()]->Write();
1090 title.str("");
1091 title << "Chi2DeltaHit_" << i + 1 << "-" << j + 1;
1092 hChi2DeltaHit[title.str()]->Write();
1093 title.str("");
1094 title << "Chi2NSharedHit_" << i + 1 << "-" << j + 1;
1095 hChi2NSharedHit[title.str()]->Write();
1096 title.str("");
1097 title << "Chi2SharedHit_" << i + 1 << "-" << j + 1;
1098 hChi2SharedHit[title.str()]->Write();
1099
1100 gp_ts->cd();
1101 gp_tsx->cd();
1102 title.str("");
1103 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_ts";
1104 hPullGP_X_ts[title.str()]->Write();
1105 gp_tsy->cd();
1106 title.str("");
1107 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_ts";
1108 hPullGP_Y_ts[title.str()]->Write();
1109 gp_tsz->cd();
1110 title.str("");
1111 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_ts";
1112 hPullGP_Z_ts[title.str()]->Write();
1113
1114 gm_ts->cd();
1115 gm_tsx->cd();
1116 title.str("");
1117 title << "PullGM_X_" << i + 1 << "-" << j + 1 << "_ts";
1118 hPullGM_X_ts[title.str()]->Write();
1119 gm_tsy->cd();
1120 title.str("");
1121 title << "PullGM_Y_" << i + 1 << "-" << j + 1 << "_ts";
1122 hPullGM_Y_ts[title.str()]->Write();
1123 gm_tsz->cd();
1124 title.str("");
1125 title << "PullGM_Z_" << i + 1 << "-" << j + 1 << "_ts";
1126 hPullGM_Z_ts[title.str()]->Write();
1127
1128 gp_tr->cd();
1129 gp_trx->cd();
1130 title.str("");
1131 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_tr";
1132 hPullGP_X_tr[title.str()]->Write();
1133 gp_try->cd();
1134 title.str("");
1135 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_tr";
1136 hPullGP_Y_tr[title.str()]->Write();
1137 gp_trz->cd();
1138 title.str("");
1139 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_tr";
1140 hPullGP_Z_tr[title.str()]->Write();
1141
1142 gp_rs->cd();
1143 gp_rsx->cd();
1144 title.str("");
1145 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_rs";
1146 hPullGP_X_rs[title.str()]->Write();
1147 gp_rsy->cd();
1148 title.str("");
1149 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_rs";
1150 hPullGP_Y_rs[title.str()]->Write();
1151 gp_rsz->cd();
1152 title.str("");
1153 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_rs";
1154 hPullGP_Z_rs[title.str()]->Write();
1155
1156 if (((i == 2 || i == 4) && (j == 0 || j == 1)) || (i == 3 || i == 5)) {
1157 chi2d->cd();
1158 title.str("");
1159 title << "Chi2Increment_mono_" << i + 1 << "-" << j + 1;
1160 hChi2Increment_mono[title.str()]->Write();
1161 title.str("");
1162 title << "Chi2Increment_stereo_" << i + 1 << "-" << j + 1;
1163 hChi2Increment_stereo[title.str()]->Write();
1164
1165 gp_ts->cd();
1166 gp_tsx_mono->cd();
1167 title.str("");
1168 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_ts_mono";
1169 hPullGP_X_ts_mono[title.str()]->Write();
1170 gp_tsy_mono->cd();
1171 title.str("");
1172 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_ts_mono";
1173 hPullGP_Y_ts_mono[title.str()]->Write();
1174 gp_tsz_mono->cd();
1175 title.str("");
1176 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_ts_mono";
1177 hPullGP_Z_ts_mono[title.str()]->Write();
1178
1179 gm_ts->cd();
1180 gm_tsx_mono->cd();
1181 title.str("");
1182 title << "PullGM_X_" << i + 1 << "-" << j + 1 << "_ts_mono";
1183 hPullGM_X_ts_mono[title.str()]->Write();
1184 gm_tsy_mono->cd();
1185 title.str("");
1186 title << "PullGM_Y_" << i + 1 << "-" << j + 1 << "_ts_mono";
1187 hPullGM_Y_ts_mono[title.str()]->Write();
1188 gm_tsz_mono->cd();
1189 title.str("");
1190 title << "PullGM_Z_" << i + 1 << "-" << j + 1 << "_ts_mono";
1191 hPullGM_Z_ts_mono[title.str()]->Write();
1192
1193 gp_tr->cd();
1194 gp_trx_mono->cd();
1195 title.str("");
1196 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_tr_mono";
1197 hPullGP_X_tr_mono[title.str()]->Write();
1198 gp_try_mono->cd();
1199 title.str("");
1200 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_tr_mono";
1201 hPullGP_Y_tr_mono[title.str()]->Write();
1202 gp_trz_mono->cd();
1203 title.str("");
1204 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_tr_mono";
1205 hPullGP_Z_tr_mono[title.str()]->Write();
1206
1207 gp_rs->cd();
1208 gp_rsx_mono->cd();
1209 title.str("");
1210 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_rs_mono";
1211 hPullGP_X_rs_mono[title.str()]->Write();
1212 gp_rsy_mono->cd();
1213 title.str("");
1214 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_rs_mono";
1215 hPullGP_Y_rs_mono[title.str()]->Write();
1216 gp_rsz_mono->cd();
1217 title.str("");
1218 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_rs_mono";
1219 hPullGP_Z_rs_mono[title.str()]->Write();
1220
1221
1222 gp_ts->cd();
1223 gp_tsx_stereo->cd();
1224 title.str("");
1225 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_ts_stereo";
1226 hPullGP_X_ts_stereo[title.str()]->Write();
1227 gp_tsy_stereo->cd();
1228 title.str("");
1229 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_ts_stereo";
1230 hPullGP_Y_ts_stereo[title.str()]->Write();
1231 gp_tsz_stereo->cd();
1232 title.str("");
1233 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_ts_stereo";
1234 hPullGP_Z_ts_stereo[title.str()]->Write();
1235
1236 gm_ts->cd();
1237 gm_tsx_stereo->cd();
1238 title.str("");
1239 title << "PullGM_X_" << i + 1 << "-" << j + 1 << "_ts_stereo";
1240 hPullGM_X_ts_stereo[title.str()]->Write();
1241 gm_tsy_stereo->cd();
1242 title.str("");
1243 title << "PullGM_Y_" << i + 1 << "-" << j + 1 << "_ts_stereo";
1244 hPullGM_Y_ts_stereo[title.str()]->Write();
1245 gm_tsz_stereo->cd();
1246 title.str("");
1247 title << "PullGM_Z_" << i + 1 << "-" << j + 1 << "_ts_stereo";
1248 hPullGM_Z_ts_stereo[title.str()]->Write();
1249
1250 gp_tr->cd();
1251 gp_trx_stereo->cd();
1252 title.str("");
1253 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_tr_stereo";
1254 hPullGP_X_tr_stereo[title.str()]->Write();
1255 gp_try_stereo->cd();
1256 title.str("");
1257 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_tr_stereo";
1258 hPullGP_Y_tr_stereo[title.str()]->Write();
1259 gp_trz_stereo->cd();
1260 title.str("");
1261 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_tr_stereo";
1262 hPullGP_Z_tr_stereo[title.str()]->Write();
1263
1264 gp_rs->cd();
1265 gp_rsx_stereo->cd();
1266 title.str("");
1267 title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_rs_stereo";
1268 hPullGP_X_rs_stereo[title.str()]->Write();
1269 gp_rsy_stereo->cd();
1270 title.str("");
1271 title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_rs_stereo";
1272 hPullGP_Y_rs_stereo[title.str()]->Write();
1273 gp_rsz_stereo->cd();
1274 title.str("");
1275 title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_rs_stereo";
1276 hPullGP_Z_rs_stereo[title.str()]->Write();
1277 }
1278 }
1279
1280 file->Close();
1281 }
1282
1283
1284
1285 std::pair<LocalPoint, LocalVector> TestTrackHits::projectHit(const PSimHit& hit,
1286 const StripGeomDetUnit* stripDet,
1287 const BoundPlane& plane) {
1288 const StripTopology& topol = stripDet->specificTopology();
1289 GlobalPoint globalpos = stripDet->surface().toGlobal(hit.localPosition());
1290 LocalPoint localHit = plane.toLocal(globalpos);
1291
1292 LocalVector locdir = hit.localDirection();
1293
1294
1295 GlobalVector globaldir = stripDet->surface().toGlobal(locdir);
1296 LocalVector dir = plane.toLocal(globaldir);
1297 float scale = -localHit.z() / dir.z();
1298
1299 LocalPoint projectedPos = localHit + scale * dir;
1300
1301 float selfAngle = topol.stripAngle(topol.strip(hit.localPosition()));
1302
1303 LocalVector stripDir(sin(selfAngle), cos(selfAngle), 0);
1304
1305 LocalVector localStripDir(plane.toLocal(stripDet->surface().toGlobal(stripDir)));
1306
1307 return std::pair<LocalPoint, LocalVector>(projectedPos, localStripDir);
1308 }
1309
1310 template <unsigned int D>
1311 double TestTrackHits::computeChi2Increment(MeasurementExtractor me, TransientTrackingRecHit::ConstRecHitPointer rhit) {
1312 typedef typename AlgebraicROOTObject<D>::Vector VecD;
1313 typedef typename AlgebraicROOTObject<D, D>::SymMatrix SMatDD;
1314 VecD r = asSVector<D>(rhit->parameters()) - me.measuredParameters<D>(*rhit);
1315
1316 SMatDD R = asSMatrix<D>(rhit->parametersError()) + me.measuredError<D>(*rhit);
1317 R.Invert();
1318 return ROOT::Math::Similarity(r, R);
1319 }
1320
1321 #include "FWCore/Framework/interface/ModuleFactory.h"
1322 #include "FWCore/Framework/interface/MakerMacros.h"
1323 DEFINE_FWK_MODULE(TestTrackHits);
