Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-01-13 23:21:01

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   //   ParameterSet cuts = iConfig.getParameter<ParameterSet>("RecoTracksCuts");
0033   //   selectRecoTracks = RecoTrackSelector(cuts.getParameter<double>("ptMin"),
0034   //                       cuts.getParameter<double>("minRapidity"),
0035   //                       cuts.getParameter<double>("maxRapidity"),
0036   //                       cuts.getParameter<double>("tip"),
0037   //                       cuts.getParameter<double>("lip"),
0038   //                       cuts.getParameter<int>("minHit"),
0039   //                       cuts.getParameter<double>("maxChi2"));
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         //mono
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         //stereo
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   //Retrieve tracker topology from geometry
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     //     if ( !selectRecoTracks( *track,beamSpot.product() ) ) {
0299     //       LogTrace("TestTrackHits") << "track does not pass quality cuts: skippingtrack #" << ++yyy;
0300     //       i++;
0301     //       continue;
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     //     if(recSimColl.find(track) != recSimColl.end()) {
0320     //       tP = recSimColl[track];
0321     //     } else {
0322     //       LogTrace("TestTrackHits") << "fake track: skipping track " << ++yyy;
0323     //       continue;//skip fake tracks
0324     //     }
0325     //     if (tP.size()==0) {
0326     //       LogTrace("TestTrackHits") << "fake track: skipping track " << ++yyy;
0327     //       continue;
0328     //     }
0329     TrackingParticleRef tp = tP.begin()->first;
0330     LogTrace("TestTrackHits") << "a tp is associated with fraction=" << tP.begin()->second;
0331     //LogTrace("TestTrackHits") << "last tp is associated with fraction=" << (tP.end()-1)->second;
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     //LogTrace("TestTrackHits") << "Analyzing hits of track number " << ++yyy << " good track number " << ++yy;
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       //TSOS state = tm->backwardPredictedState();
0346       //TSOS state = tm->forwardPredictedState();
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           //std::cout << "detUnitId="  << m->detUnitId() << " trackId=" << m->trackId() << " energyLoss=" << m->energyLoss() << std::endl;
0392         } else {
0393           energyLossM.push_back(m->energyLoss());
0394         }
0395       }
0396       //double delta = 99999;
0397       //LocalPoint rhitLPv = rhit->localPosition();
0398       //vector<PSimHit> assSimHits = hitAssociator.associateHit(*(rhit)->hit());
0399       //unsigned int  simhitvecsize = assSimHits.size();
0400       //if (simhitvecsize==0) continue;
0401       //PSimHit shit;
0402       //for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
0403       //  if ((m->localPosition()-rhitLPv).mag()<delta) {
0404       //    shit=*m;
0405       //    delta = (m->localPosition()-rhitLPv).mag();
0406       //  }
0407       //}
0408 
0409       //plot chi2 increment
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       //       if (dynamic_cast<const SiPixelRecHit*>(rhit->hit()))
0462       //    hClsize_vs_Chi2->Fill( chi2increment, ((const SiPixelRecHit*)(rhit->hit()))->cluster()->size() );
0463       //       if (dynamic_cast<const SiStripRecHit2D*>(rhit->hit()))
0464       //    hClsize_vs_Chi2->Fill( chi2increment, ((const SiStripRecHit2D*)(rhit->hit()))->cluster()->amplitudes().size() );
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           //not optimized for matched hits
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;  //ger.similarity(delta);
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           //project simhit;
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       //double pullGPX_rs = (rhitGPos.x()-shitGPos.x());
0678       //double pullGPY_rs = (rhitGPos.y()-shitGPos.y());
0679       //double pullGPZ_rs = (rhitGPos.z()-shitGPos.z());
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       //double pullGPX_tr = (tsosGPos.x()-rhitGPos.x());
0697       //double pullGPY_tr = (tsosGPos.y()-rhitGPos.y());
0698       //double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z());
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       //double pullGPX_ts = (tsosGPos.x()-shitGPos.x());
0716       //double pullGPY_ts = (tsosGPos.y()-shitGPos.y());
0717       //double pullGPZ_ts = (tsosGPos.z()-shitGPos.z());
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       //double pullGMX_ts = (tsosGMom.x()-shitGMom.x());
0735       //double pullGMY_ts = (tsosGMom.y()-shitGMom.y());
0736       //double pullGMZ_ts = (tsosGMom.z()-shitGMom.z());
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         //mono
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         //LogTrace("TestTrackHits") << "assMonoSimHits.size()=" << assMonoSimHits.size() ;
0774         //LogTrace("TestTrackHits") << "mono shit=" << monoShitGPos ;
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         //double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x());
0788         //double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y());
0789         //double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z());
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         //double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x());
0812         //double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y());
0813         //double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z());
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         //double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x());
0829         //double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y());
0830         //double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z());
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         //double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x());
0846         //double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y());
0847         //double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z());
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         //stereo
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         //LogTrace("TestTrackHits") << "assStereoSimHits.size()=" << assStereoSimHits.size() ;
0882         //LogTrace("TestTrackHits") << "stereo shit=" << stereoShitGPos ;
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         //  double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x())/*/sqrt(stereoRhitGPEr.cxx())*/;
0903         //  double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y())/*/sqrt(stereoRhitGPEr.cyy())*/;
0904         //  double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z())/*/sqrt(stereoRhitGPEr.czz())*/;
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         //double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x());
0923         //double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y());
0924         //double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z());
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         //double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x());
0940         //double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y());
0941         //double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z());
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         //double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x());
0957         //double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y());
0958         //double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z());
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   //file->Write();
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         //mono
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         //stereo
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 //needed by to do the residual for matched hits
1284 //taken from SiStripTrackingRecHitsValid.cc
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   //track direction
1292   LocalVector locdir = hit.localDirection();
1293   //rotate track in new frame
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);  // vector along strip in hit frame
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);