Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #include "TestSmoothHits.h"
0002 
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0005 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0006 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0007 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0008 #include "DataFormats/TrackerRecHit2D/interface/TrackingRecHitLessFromGlobalPosition.h"
0009 #include "TrackingTools/PatternTools/interface/TSCPBuilderNoMaterial.h"
0010 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0011 #include <TDirectory.h>
0012 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0013 
0014 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0015 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0016 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0017 
0018 typedef TrajectoryStateOnSurface TSOS;
0019 typedef TransientTrackingRecHit::ConstRecHitPointer CTTRHp;
0020 using namespace std;
0021 using namespace edm;
0022 
0023 TestSmoothHits::TestSmoothHits(const edm::ParameterSet& iConfig) : trackerHitAssociatorConfig_(consumesCollector()) {
0024   LogTrace("TestSmoothHits") << iConfig << std::endl;
0025   propagatorName = iConfig.getParameter<std::string>("Propagator");
0026   builderName = iConfig.getParameter<std::string>("TTRHBuilder");
0027   srcName = iConfig.getParameter<std::string>("src");
0028   fname = iConfig.getParameter<std::string>("Fitter");
0029   sname = iConfig.getParameter<std::string>("Smoother");
0030   mineta = iConfig.getParameter<double>("mineta");
0031   maxeta = iConfig.getParameter<double>("maxeta");
0032 
0033   theGToken = esConsumes<edm::Transition::BeginRun>();
0034   theMFToken = esConsumes<edm::Transition::BeginRun>();
0035   thePropagatorToken = esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", propagatorName));
0036   theBuilderToken = esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", builderName));
0037   fitToken = esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", fname));
0038   smoothToken = esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", sname));
0039   tTopoToken = esConsumes();
0040 
0041   theTCCollectionToken = consumes(edm::InputTag(srcName));
0042 }
0043 
0044 TestSmoothHits::~TestSmoothHits() {}
0045 
0046 void TestSmoothHits::beginRun(edm::Run const& run, const edm::EventSetup& iSetup) {
0047   theG = iSetup.getHandle(theGToken);
0048   theMF = iSetup.getHandle(theMFToken);
0049   thePropagator = iSetup.getHandle(thePropagatorToken);
0050   theBuilder = iSetup.getHandle(theBuilderToken);
0051   fit = iSetup.getHandle(fitToken);
0052   smooth = iSetup.getHandle(smoothToken);
0053 
0054   file = new TFile("testSmoothHits.root", "recreate");
0055   for (int i = 0; i != 6; i++)
0056     for (int j = 0; j != 9; j++) {
0057       if (i == 0 && j > 2)
0058         break;
0059       if (i == 1 && j > 1)
0060         break;
0061       if (i == 2 && j > 3)
0062         break;
0063       if (i == 3 && j > 2)
0064         break;
0065       if (i == 4 && j > 5)
0066         break;
0067       if (i == 5 && j > 8)
0068         break;
0069       title.str("");
0070       title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_ts";
0071       hPullGP_X_ts[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0072       title.str("");
0073       title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_ts";
0074       hPullGP_Y_ts[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0075       title.str("");
0076       title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_ts";
0077       hPullGP_Z_ts[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0078       title.str("");
0079       title << "Chi2Increment_" << i + 1 << "-" << j + 1;
0080       hChi2Increment[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, 0, 100);
0081 
0082       title.str("");
0083       title << "PullGM_X_" << i + 1 << "-" << j + 1 << "_ts";
0084       hPullGM_X_ts[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0085       title.str("");
0086       title << "PullGM_Y_" << i + 1 << "-" << j + 1 << "_ts";
0087       hPullGM_Y_ts[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0088       title.str("");
0089       title << "PullGM_Z_" << i + 1 << "-" << j + 1 << "_ts";
0090       hPullGM_Z_ts[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0091 
0092       title.str("");
0093       title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_tr";
0094       hPullGP_X_tr[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0095       title.str("");
0096       title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_tr";
0097       hPullGP_Y_tr[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0098       title.str("");
0099       title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_tr";
0100       hPullGP_Z_tr[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0101 
0102       title.str("");
0103       title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_rs";
0104       hPullGP_X_rs[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0105       title.str("");
0106       title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_rs";
0107       hPullGP_Y_rs[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0108       title.str("");
0109       title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_rs";
0110       hPullGP_Z_rs[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0111 
0112       if (((i == 2 || i == 4) && (j == 0 || j == 1)) || (i == 3 || i == 5)) {
0113         //mono
0114         title.str("");
0115         title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_ts_mono";
0116         hPullGP_X_ts_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0117         title.str("");
0118         title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_ts_mono";
0119         hPullGP_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0120         title.str("");
0121         title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_ts_mono";
0122         hPullGP_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0123 
0124         title.str("");
0125         title << "PullGM_X_" << i + 1 << "-" << j + 1 << "_ts_mono";
0126         hPullGM_X_ts_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0127         title.str("");
0128         title << "PullGM_Y_" << i + 1 << "-" << j + 1 << "_ts_mono";
0129         hPullGM_Y_ts_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0130         title.str("");
0131         title << "PullGM_Z_" << i + 1 << "-" << j + 1 << "_ts_mono";
0132         hPullGM_Z_ts_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0133 
0134         title.str("");
0135         title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_tr_mono";
0136         hPullGP_X_tr_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0137         title.str("");
0138         title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_tr_mono";
0139         hPullGP_Y_tr_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0140         title.str("");
0141         title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_tr_mono";
0142         hPullGP_Z_tr_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0143 
0144         title.str("");
0145         title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_rs_mono";
0146         hPullGP_X_rs_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0147         title.str("");
0148         title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_rs_mono";
0149         hPullGP_Y_rs_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0150         title.str("");
0151         title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_rs_mono";
0152         hPullGP_Z_rs_mono[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0153 
0154         //stereo
0155         title.str("");
0156         title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0157         hPullGP_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0158         title.str("");
0159         title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0160         hPullGP_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0161         title.str("");
0162         title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0163         hPullGP_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0164 
0165         title.str("");
0166         title << "PullGM_X_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0167         hPullGM_X_ts_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0168         title.str("");
0169         title << "PullGM_Y_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0170         hPullGM_Y_ts_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0171         title.str("");
0172         title << "PullGM_Z_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0173         hPullGM_Z_ts_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0174 
0175         title.str("");
0176         title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_tr_stereo";
0177         hPullGP_X_tr_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0178         title.str("");
0179         title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_tr_stereo";
0180         hPullGP_Y_tr_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0181         title.str("");
0182         title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_tr_stereo";
0183         hPullGP_Z_tr_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0184 
0185         title.str("");
0186         title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_rs_stereo";
0187         hPullGP_X_rs_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0188         title.str("");
0189         title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_rs_stereo";
0190         hPullGP_Y_rs_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0191         title.str("");
0192         title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_rs_stereo";
0193         hPullGP_Z_rs_stereo[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0194       }
0195     }
0196   hTotChi2Increment = new TH1F("TotChi2Increment", "TotChi2Increment", 1000, 0, 100);
0197   hChi2_vs_Process = new TH2F("Chi2_vs_Process", "Chi2_vs_Process", 1000, 0, 100, 17, -0.5, 16.5);
0198   hChi2_vs_clsize = new TH2F("Chi2_vs_clsize", "Chi2_vs_clsize", 1000, 0, 100, 17, -0.5, 16.5);
0199 }
0200 
0201 void TestSmoothHits::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0202   //Retrieve tracker topology from geometry
0203   edm::ESHandle<TrackerTopology> tTopo = iSetup.getHandle(tTopoToken);
0204 
0205   LogTrace("TestSmoothHits") << "new event" << std::endl;
0206 
0207   theTCCollection = iEvent.getHandle(theTCCollectionToken);
0208   TrackerHitAssociator hitAssociator(iEvent, trackerHitAssociatorConfig_);
0209 
0210   TrajectoryStateCombiner combiner;
0211 
0212   for (TrackCandidateCollection::const_iterator i = theTCCollection->begin(); i != theTCCollection->end(); i++) {
0213     LogTrace("TestSmoothHits") << "new candidate" << std::endl;
0214 
0215     const TrackCandidate* theTC = &(*i);
0216     PTrajectoryStateOnDet state = theTC->trajectoryStateOnDet();
0217 
0218     //convert PTrajectoryStateOnDet to TrajectoryStateOnSurface
0219 
0220     DetId detId(state.detId());
0221     TrajectoryStateOnSurface theTSOS =
0222         trajectoryStateTransform::transientState(state, &(theG->idToDet(detId)->surface()), theMF.product());
0223 
0224     if (theTSOS.globalMomentum().eta() > maxeta || theTSOS.globalMomentum().eta() < mineta)
0225       continue;
0226 
0227     //convert the TrackingRecHit vector to a TransientTrackingRecHit vector
0228     TransientTrackingRecHit::RecHitContainer hits;
0229 
0230     for (auto const& recHit : theTC->recHits()) {
0231       hits.push_back(theBuilder->build(&recHit));
0232     }
0233 
0234     //call the fitter
0235     std::vector<Trajectory> fitted = fit->fit(theTC->seed(), hits, theTSOS);
0236     //call the smoother
0237     std::vector<Trajectory> result;
0238     for (std::vector<Trajectory>::iterator it = fitted.begin(); it != fitted.end(); it++) {
0239       std::vector<Trajectory> smoothed = smooth->trajectories(*it);
0240       result.insert(result.end(), smoothed.begin(), smoothed.end());
0241     }
0242     if (result.empty())
0243       continue;
0244     std::vector<TrajectoryMeasurement> vtm = result[0].measurements();
0245 
0246     TSOS lastState = theTSOS;
0247     for (std::vector<TrajectoryMeasurement>::iterator tm = vtm.begin(); tm != vtm.end(); tm++) {
0248       TransientTrackingRecHit::ConstRecHitPointer rhit = tm->recHit();
0249       if ((rhit)->isValid() == 0 && rhit->det() != nullptr)
0250         continue;
0251       LogTrace("TestSmoothHits") << "new hit";
0252 
0253       int subdetId = rhit->det()->geographicalId().subdetId();
0254       DetId id = rhit->det()->geographicalId();
0255       int layerId = tTopo->layer(id);
0256       LogTrace("TestSmoothHits") << "subdetId=" << subdetId << " layerId=" << layerId;
0257 
0258       double delta = 99999;
0259       LocalPoint rhitLPv = rhit->localPosition();
0260 
0261       std::vector<PSimHit> assSimHits = hitAssociator.associateHit(*(rhit->hit()));
0262       if (assSimHits.empty())
0263         continue;
0264       PSimHit shit;
0265       for (std::vector<PSimHit>::const_iterator m = assSimHits.begin(); m < assSimHits.end(); m++) {
0266         if ((m->localPosition() - rhitLPv).mag() < delta) {
0267           shit = *m;
0268           delta = (m->localPosition() - rhitLPv).mag();
0269         }
0270       }
0271 
0272       TSOS currentState = combiner(tm->backwardPredictedState(), tm->forwardPredictedState());
0273       TSOS updatedState = tm->updatedState();
0274 
0275       //plot chi2 increment
0276       double chi2increment = tm->estimate();
0277       LogTrace("TestSmoothHits") << "tm->estimate()=" << tm->estimate();
0278       title.str("");
0279       title << "Chi2Increment_" << subdetId << "-" << layerId;
0280       hChi2Increment[title.str()]->Fill(chi2increment);
0281       hTotChi2Increment->Fill(chi2increment);
0282       hChi2_vs_Process->Fill(chi2increment, shit.processType());
0283       if (dynamic_cast<const SiPixelRecHit*>(rhit->hit()))
0284         hChi2_vs_clsize->Fill(chi2increment, ((const SiPixelRecHit*)(rhit->hit()))->cluster()->size());
0285       if (dynamic_cast<const SiStripRecHit2D*>(rhit->hit()))
0286         hChi2_vs_clsize->Fill(chi2increment, ((const SiStripRecHit2D*)(rhit->hit()))->cluster()->amplitudes().size());
0287 
0288       //test hits
0289       const Surface* surf = &((rhit)->det()->surface());
0290       LocalVector shitLMom;
0291       LocalPoint shitLPos;
0292       if (dynamic_cast<const SiStripMatchedRecHit2D*>(rhit->hit())) {
0293         double rechitmatchedx = rhit->localPosition().x();
0294         double rechitmatchedy = rhit->localPosition().y();
0295         double mindist = 999999;
0296         float distx, disty;
0297         std::pair<LocalPoint, LocalVector> closestPair;
0298         const StripGeomDetUnit* stripDet = (StripGeomDetUnit*)((const GluedGeomDet*)(rhit)->det())->stereoDet();
0299         const BoundPlane& plane = (rhit)->det()->surface();
0300         for (std::vector<PSimHit>::const_iterator m = assSimHits.begin(); m < assSimHits.end(); m++) {
0301           //project simhit;
0302           std::pair<LocalPoint, LocalVector> hitPair = projectHit((*m), stripDet, plane);
0303           distx = fabs(rechitmatchedx - hitPair.first.x());
0304           disty = fabs(rechitmatchedy - hitPair.first.y());
0305           double dist = distx * distx + disty * disty;
0306           if (sqrt(dist) < mindist) {
0307             mindist = dist;
0308             closestPair = hitPair;
0309           }
0310         }
0311         shitLPos = closestPair.first;
0312         shitLMom = closestPair.second;
0313       } else {
0314         shitLPos = shit.localPosition();
0315         shitLMom = shit.momentumAtEntry();
0316       }
0317       GlobalVector shitGMom = surf->toGlobal(shitLMom);
0318       GlobalPoint shitGPos = surf->toGlobal(shitLPos);
0319       //      if (dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())) {
0320       //    double rechitmatchedx = rhit->localPosition().x();
0321       //    double rechitmatchedy = rhit->localPosition().y();
0322       //    double mindist = 999999;
0323       //    double distx, disty;
0324       //    std::pair<LocalPoint,LocalVector> closestPair;
0325       //    const StripGeomDetUnit* stripDet =(StripGeomDetUnit*) ((const GluedGeomDet *)(rhit)->det())->stereoDet();
0326       //    const BoundPlane& plane = (rhit)->det()->surface();
0327       //    for(std::vector<PSimHit>::const_iterator m=assSimHits.begin(); m<assSimHits.end(); m++){
0328       //      const PSimHit& hit = *m;
0329       //      const StripTopology& topol = stripDet->specificTopology();
0330       //      GlobalPoint globalpos = stripDet->surface().toGlobal(hit.localPosition());
0331       //      LocalPoint localHit = plane.toLocal(globalpos);
0332       //      //track direction
0333       //      LocalVector locdir=hit.localDirection();
0334       //      //rotate track in new frame
0335       //      GlobalVector globaldir= stripDet->surface().toGlobal(locdir);
0336       //      LocalVector dir=plane.toLocal(globaldir);
0337       //      float scale = -localHit.z() / dir.z();
0338       //      LocalPoint projectedPos = localHit + scale*dir;
0339       //      float selfAngle = topol.stripAngle( topol.strip( hit.localPosition()));
0340       //      LocalVector stripDir( sin(selfAngle), cos(selfAngle), 0); // vector along strip in hit frame
0341       //      LocalVector localStripDir = LocalVector( plane.toLocal(stripDet->surface().toGlobal( stripDir)));
0342       //      std::pair<LocalPoint,LocalVector> hitPair( projectedPos, localStripDir);
0343       //      distx = fabs(rechitmatchedx - hitPair.first.x());
0344       //      disty = fabs(rechitmatchedy - hitPair.first.y());
0345       //      double dist = distx*distx+disty*disty;
0346       //      if(sqrt(dist)<mindist){
0347       //        mindist = dist;
0348       //        closestPair = hitPair;
0349       //      }
0350       //    }
0351       //    shitLPos = closestPair.first;
0352       //    shitLMom = closestPair.second;
0353       //       } else {
0354       //    shitLPos = shit.localPosition();
0355       //    shitLMom = shit.momentumAtEntry();
0356       //       }
0357       //       GlobalVector shitGMom = surf->toGlobal(shitLMom);
0358       //       GlobalPoint shitGPos = surf->toGlobal(shitLPos);
0359 
0360       GlobalVector tsosGMom = currentState.globalMomentum();
0361       GlobalError tsosGMEr(currentState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3, 3));
0362       GlobalPoint tsosGPos = currentState.globalPosition();
0363       GlobalError tsosGPEr = currentState.cartesianError().position();
0364 
0365       GlobalPoint rhitGPos = (rhit)->globalPosition();
0366       GlobalError rhitGPEr = (rhit)->globalPositionError();
0367 
0368       double pullGPX_rs = (rhitGPos.x() - shitGPos.x()) / sqrt(rhitGPEr.cxx());
0369       double pullGPY_rs = (rhitGPos.y() - shitGPos.y()) / sqrt(rhitGPEr.cyy());
0370       double pullGPZ_rs = (rhitGPos.z() - shitGPos.z()) / sqrt(rhitGPEr.czz());
0371       //double pullGPX_rs = (rhitGPos.x()-shitGPos.x());
0372       //double pullGPY_rs = (rhitGPos.y()-shitGPos.y());
0373       //double pullGPZ_rs = (rhitGPos.z()-shitGPos.z());
0374 
0375       LogTrace("TestSmoothHits") << "rs" << std::endl;
0376 
0377       title.str("");
0378       title << "PullGP_X_" << subdetId << "-" << layerId << "_rs";
0379       hPullGP_X_rs[title.str()]->Fill(pullGPX_rs);
0380       title.str("");
0381       title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs";
0382       hPullGP_Y_rs[title.str()]->Fill(pullGPY_rs);
0383       title.str("");
0384       title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs";
0385       hPullGP_Z_rs[title.str()]->Fill(pullGPZ_rs);
0386 
0387       double pullGPX_tr = (tsosGPos.x() - rhitGPos.x()) / sqrt(tsosGPEr.cxx() + rhitGPEr.cxx());
0388       double pullGPY_tr = (tsosGPos.y() - rhitGPos.y()) / sqrt(tsosGPEr.cyy() + rhitGPEr.cyy());
0389       double pullGPZ_tr = (tsosGPos.z() - rhitGPos.z()) / sqrt(tsosGPEr.czz() + rhitGPEr.czz());
0390       //double pullGPX_tr = (tsosGPos.x()-rhitGPos.x());
0391       //double pullGPY_tr = (tsosGPos.y()-rhitGPos.y());
0392       //double pullGPZ_tr = (tsosGPos.z()-rhitGPos.z());
0393 
0394       LogTrace("TestSmoothHits") << "tr" << std::endl;
0395 
0396       title.str("");
0397       title << "PullGP_X_" << subdetId << "-" << layerId << "_tr";
0398       hPullGP_X_tr[title.str()]->Fill(pullGPX_tr);
0399       title.str("");
0400       title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr";
0401       hPullGP_Y_tr[title.str()]->Fill(pullGPY_tr);
0402       title.str("");
0403       title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr";
0404       hPullGP_Z_tr[title.str()]->Fill(pullGPZ_tr);
0405 
0406       double pullGPX_ts = (tsosGPos.x() - shitGPos.x()) / sqrt(tsosGPEr.cxx());
0407       double pullGPY_ts = (tsosGPos.y() - shitGPos.y()) / sqrt(tsosGPEr.cyy());
0408       double pullGPZ_ts = (tsosGPos.z() - shitGPos.z()) / sqrt(tsosGPEr.czz());
0409       //double pullGPX_ts = (tsosGPos.x()-shitGPos.x());
0410       //double pullGPY_ts = (tsosGPos.y()-shitGPos.y());
0411       //double pullGPZ_ts = (tsosGPos.z()-shitGPos.z());
0412 
0413       LogTrace("TestSmoothHits") << "ts1" << std::endl;
0414 
0415       title.str("");
0416       title << "PullGP_X_" << subdetId << "-" << layerId << "_ts";
0417       hPullGP_X_ts[title.str()]->Fill(pullGPX_ts);
0418       title.str("");
0419       title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts";
0420       hPullGP_Y_ts[title.str()]->Fill(pullGPY_ts);
0421       title.str("");
0422       title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts";
0423       hPullGP_Z_ts[title.str()]->Fill(pullGPZ_ts);
0424 
0425       double pullGMX_ts = (tsosGMom.x() - shitGMom.x()) / sqrt(tsosGMEr.cxx());
0426       double pullGMY_ts = (tsosGMom.y() - shitGMom.y()) / sqrt(tsosGMEr.cyy());
0427       double pullGMZ_ts = (tsosGMom.z() - shitGMom.z()) / sqrt(tsosGMEr.czz());
0428       //double pullGMX_ts = (tsosGMom.x()-shitGMom.x());
0429       //double pullGMY_ts = (tsosGMom.y()-shitGMom.y());
0430       //double pullGMZ_ts = (tsosGMom.z()-shitGMom.z());
0431 
0432       LogTrace("TestSmoothHits") << "ts2" << std::endl;
0433 
0434       title.str("");
0435       title << "PullGM_X_" << subdetId << "-" << layerId << "_ts";
0436       hPullGM_X_ts[title.str()]->Fill(pullGMX_ts);
0437       title.str("");
0438       title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts";
0439       hPullGM_Y_ts[title.str()]->Fill(pullGMY_ts);
0440       title.str("");
0441       title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts";
0442       hPullGM_Z_ts[title.str()]->Fill(pullGMZ_ts);
0443 
0444       if (dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())) {
0445         //mono
0446         LogTrace("TestSmoothHits") << "MONO HIT" << std::endl;
0447         auto m = dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())->monoHit();
0448         CTTRHp tMonoHit = theBuilder->build(&m);
0449         if (tMonoHit == nullptr)
0450           continue;
0451         vector<PSimHit> assMonoSimHits = hitAssociator.associateHit(*tMonoHit->hit());
0452         if (assMonoSimHits.empty())
0453           continue;
0454         const PSimHit sMonoHit = *(assSimHits.begin());
0455         const Surface* monoSurf = &(tMonoHit->det()->surface());
0456         if (monoSurf == nullptr)
0457           continue;
0458         TSOS monoState = thePropagator->propagate(lastState, *monoSurf);
0459         if (monoState.isValid() == 0)
0460           continue;
0461 
0462         LocalVector monoShitLMom = sMonoHit.momentumAtEntry();
0463         GlobalVector monoShitGMom = monoSurf->toGlobal(monoShitLMom);
0464         LocalPoint monoShitLPos = sMonoHit.localPosition();
0465         GlobalPoint monoShitGPos = monoSurf->toGlobal(monoShitLPos);
0466 
0467         GlobalVector monoTsosGMom = monoState.globalMomentum();
0468         GlobalError monoTsosGMEr(monoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3, 3));
0469         GlobalPoint monoTsosGPos = monoState.globalPosition();
0470         GlobalError monoTsosGPEr = monoState.cartesianError().position();
0471 
0472         GlobalPoint monoRhitGPos = tMonoHit->globalPosition();
0473         GlobalError monoRhitGPEr = tMonoHit->globalPositionError();
0474 
0475         double pullGPX_rs_mono = (monoRhitGPos.x() - monoShitGPos.x()) / sqrt(monoRhitGPEr.cxx());
0476         double pullGPY_rs_mono = (monoRhitGPos.y() - monoShitGPos.y()) / sqrt(monoRhitGPEr.cyy());
0477         double pullGPZ_rs_mono = (monoRhitGPos.z() - monoShitGPos.z()) / sqrt(monoRhitGPEr.czz());
0478         //double pullGPX_rs_mono = (monoRhitGPos.x()-monoShitGPos.x());
0479         //double pullGPY_rs_mono = (monoRhitGPos.y()-monoShitGPos.y());
0480         //double pullGPZ_rs_mono = (monoRhitGPos.z()-monoShitGPos.z());
0481 
0482         title.str("");
0483         title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_mono";
0484         hPullGP_X_rs_mono[title.str()]->Fill(pullGPX_rs_mono);
0485         title.str("");
0486         title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_mono";
0487         hPullGP_Y_rs_mono[title.str()]->Fill(pullGPY_rs_mono);
0488         title.str("");
0489         title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_mono";
0490         hPullGP_Z_rs_mono[title.str()]->Fill(pullGPZ_rs_mono);
0491 
0492         double pullGPX_tr_mono = (monoTsosGPos.x() - monoRhitGPos.x()) / sqrt(monoTsosGPEr.cxx() + monoRhitGPEr.cxx());
0493         double pullGPY_tr_mono = (monoTsosGPos.y() - monoRhitGPos.y()) / sqrt(monoTsosGPEr.cyy() + monoRhitGPEr.cyy());
0494         double pullGPZ_tr_mono = (monoTsosGPos.z() - monoRhitGPos.z()) / sqrt(monoTsosGPEr.czz() + monoRhitGPEr.czz());
0495         //double pullGPX_tr_mono = (monoTsosGPos.x()-monoRhitGPos.x());
0496         //double pullGPY_tr_mono = (monoTsosGPos.y()-monoRhitGPos.y());
0497         //double pullGPZ_tr_mono = (monoTsosGPos.z()-monoRhitGPos.z());
0498 
0499         title.str("");
0500         title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_mono";
0501         hPullGP_X_tr_mono[title.str()]->Fill(pullGPX_tr_mono);
0502         title.str("");
0503         title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_mono";
0504         hPullGP_Y_tr_mono[title.str()]->Fill(pullGPY_tr_mono);
0505         title.str("");
0506         title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_mono";
0507         hPullGP_Z_tr_mono[title.str()]->Fill(pullGPZ_tr_mono);
0508 
0509         double pullGPX_ts_mono = (monoTsosGPos.x() - monoShitGPos.x()) / sqrt(monoTsosGPEr.cxx());
0510         double pullGPY_ts_mono = (monoTsosGPos.y() - monoShitGPos.y()) / sqrt(monoTsosGPEr.cyy());
0511         double pullGPZ_ts_mono = (monoTsosGPos.z() - monoShitGPos.z()) / sqrt(monoTsosGPEr.czz());
0512         //double pullGPX_ts_mono = (monoTsosGPos.x()-monoShitGPos.x());
0513         //double pullGPY_ts_mono = (monoTsosGPos.y()-monoShitGPos.y());
0514         //double pullGPZ_ts_mono = (monoTsosGPos.z()-monoShitGPos.z());
0515 
0516         title.str("");
0517         title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_mono";
0518         hPullGP_X_ts_mono[title.str()]->Fill(pullGPX_ts_mono);
0519         title.str("");
0520         title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_mono";
0521         hPullGP_Y_ts_mono[title.str()]->Fill(pullGPY_ts_mono);
0522         title.str("");
0523         title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_mono";
0524         hPullGP_Z_ts_mono[title.str()]->Fill(pullGPZ_ts_mono);
0525 
0526         double pullGMX_ts_mono = (monoTsosGMom.x() - monoShitGMom.x()) / sqrt(monoTsosGMEr.cxx());
0527         double pullGMY_ts_mono = (monoTsosGMom.y() - monoShitGMom.y()) / sqrt(monoTsosGMEr.cyy());
0528         double pullGMZ_ts_mono = (monoTsosGMom.z() - monoShitGMom.z()) / sqrt(monoTsosGMEr.czz());
0529         //double pullGMX_ts_mono = (monoTsosGMom.x()-monoShitGMom.x());
0530         //double pullGMY_ts_mono = (monoTsosGMom.y()-monoShitGMom.y());
0531         //double pullGMZ_ts_mono = (monoTsosGMom.z()-monoShitGMom.z());
0532 
0533         title.str("");
0534         title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_mono";
0535         hPullGM_X_ts_mono[title.str()]->Fill(pullGMX_ts_mono);
0536         title.str("");
0537         title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_mono";
0538         hPullGM_Y_ts_mono[title.str()]->Fill(pullGMY_ts_mono);
0539         title.str("");
0540         title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_mono";
0541         hPullGM_Z_ts_mono[title.str()]->Fill(pullGMZ_ts_mono);
0542 
0543         //stereo
0544         LogTrace("TestSmoothHits") << "STEREO HIT" << std::endl;
0545         auto s = dynamic_cast<const SiStripMatchedRecHit2D*>((rhit)->hit())->stereoHit();
0546         CTTRHp tStereoHit = theBuilder->build(&s);
0547         if (tStereoHit == nullptr)
0548           continue;
0549         vector<PSimHit> assStereoSimHits = hitAssociator.associateHit(*tStereoHit->hit());
0550         if (assStereoSimHits.empty())
0551           continue;
0552         const PSimHit sStereoHit = *(assSimHits.begin());
0553         const Surface* stereoSurf = &(tStereoHit->det()->surface());
0554         if (stereoSurf == nullptr)
0555           continue;
0556         TSOS stereoState = thePropagator->propagate(lastState, *stereoSurf);
0557         if (stereoState.isValid() == 0)
0558           continue;
0559 
0560         LocalVector stereoShitLMom = sStereoHit.momentumAtEntry();
0561         GlobalVector stereoShitGMom = stereoSurf->toGlobal(stereoShitLMom);
0562         LocalPoint stereoShitLPos = sStereoHit.localPosition();
0563         GlobalPoint stereoShitGPos = stereoSurf->toGlobal(stereoShitLPos);
0564 
0565         GlobalVector stereoTsosGMom = stereoState.globalMomentum();
0566         GlobalError stereoTsosGMEr(stereoState.cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3, 3));
0567         GlobalPoint stereoTsosGPos = stereoState.globalPosition();
0568         GlobalError stereoTsosGPEr = stereoState.cartesianError().position();
0569 
0570         GlobalPoint stereoRhitGPos = tStereoHit->globalPosition();
0571         GlobalError stereoRhitGPEr = tStereoHit->globalPositionError();
0572 
0573         double pullGPX_rs_stereo = (stereoRhitGPos.x() - stereoShitGPos.x()) / sqrt(stereoRhitGPEr.cxx());
0574         double pullGPY_rs_stereo = (stereoRhitGPos.y() - stereoShitGPos.y()) / sqrt(stereoRhitGPEr.cyy());
0575         double pullGPZ_rs_stereo = (stereoRhitGPos.z() - stereoShitGPos.z()) / sqrt(stereoRhitGPEr.czz());
0576         //double pullGPX_rs_stereo = (stereoRhitGPos.x()-stereoShitGPos.x());
0577         //double pullGPY_rs_stereo = (stereoRhitGPos.y()-stereoShitGPos.y());
0578         //double pullGPZ_rs_stereo = (stereoRhitGPos.z()-stereoShitGPos.z());
0579 
0580         title.str("");
0581         title << "PullGP_X_" << subdetId << "-" << layerId << "_rs_stereo";
0582         hPullGP_X_rs_stereo[title.str()]->Fill(pullGPX_rs_stereo);
0583         title.str("");
0584         title << "PullGP_Y_" << subdetId << "-" << layerId << "_rs_stereo";
0585         hPullGP_Y_rs_stereo[title.str()]->Fill(pullGPY_rs_stereo);
0586         title.str("");
0587         title << "PullGP_Z_" << subdetId << "-" << layerId << "_rs_stereo";
0588         hPullGP_Z_rs_stereo[title.str()]->Fill(pullGPZ_rs_stereo);
0589 
0590         double pullGPX_tr_stereo =
0591             (stereoTsosGPos.x() - stereoRhitGPos.x()) / sqrt(stereoTsosGPEr.cxx() + stereoRhitGPEr.cxx());
0592         double pullGPY_tr_stereo =
0593             (stereoTsosGPos.y() - stereoRhitGPos.y()) / sqrt(stereoTsosGPEr.cyy() + stereoRhitGPEr.cyy());
0594         double pullGPZ_tr_stereo =
0595             (stereoTsosGPos.z() - stereoRhitGPos.z()) / sqrt(stereoTsosGPEr.czz() + stereoRhitGPEr.czz());
0596         //double pullGPX_tr_stereo = (stereoTsosGPos.x()-stereoRhitGPos.x());
0597         //double pullGPY_tr_stereo = (stereoTsosGPos.y()-stereoRhitGPos.y());
0598         //double pullGPZ_tr_stereo = (stereoTsosGPos.z()-stereoRhitGPos.z());
0599 
0600         title.str("");
0601         title << "PullGP_X_" << subdetId << "-" << layerId << "_tr_stereo";
0602         hPullGP_X_tr_stereo[title.str()]->Fill(pullGPX_tr_stereo);
0603         title.str("");
0604         title << "PullGP_Y_" << subdetId << "-" << layerId << "_tr_stereo";
0605         hPullGP_Y_tr_stereo[title.str()]->Fill(pullGPY_tr_stereo);
0606         title.str("");
0607         title << "PullGP_Z_" << subdetId << "-" << layerId << "_tr_stereo";
0608         hPullGP_Z_tr_stereo[title.str()]->Fill(pullGPZ_tr_stereo);
0609 
0610         double pullGPX_ts_stereo = (stereoTsosGPos.x() - stereoShitGPos.x()) / sqrt(stereoTsosGPEr.cxx());
0611         double pullGPY_ts_stereo = (stereoTsosGPos.y() - stereoShitGPos.y()) / sqrt(stereoTsosGPEr.cyy());
0612         double pullGPZ_ts_stereo = (stereoTsosGPos.z() - stereoShitGPos.z()) / sqrt(stereoTsosGPEr.czz());
0613         //double pullGPX_ts_stereo = (stereoTsosGPos.x()-stereoShitGPos.x());
0614         //double pullGPY_ts_stereo = (stereoTsosGPos.y()-stereoShitGPos.y());
0615         //double pullGPZ_ts_stereo = (stereoTsosGPos.z()-stereoShitGPos.z());
0616 
0617         title.str("");
0618         title << "PullGP_X_" << subdetId << "-" << layerId << "_ts_stereo";
0619         hPullGP_X_ts_stereo[title.str()]->Fill(pullGPX_ts_stereo);
0620         title.str("");
0621         title << "PullGP_Y_" << subdetId << "-" << layerId << "_ts_stereo";
0622         hPullGP_Y_ts_stereo[title.str()]->Fill(pullGPY_ts_stereo);
0623         title.str("");
0624         title << "PullGP_Z_" << subdetId << "-" << layerId << "_ts_stereo";
0625         hPullGP_Z_ts_stereo[title.str()]->Fill(pullGPZ_ts_stereo);
0626 
0627         double pullGMX_ts_stereo = (stereoTsosGMom.x() - stereoShitGMom.x()) / sqrt(stereoTsosGMEr.cxx());
0628         double pullGMY_ts_stereo = (stereoTsosGMom.y() - stereoShitGMom.y()) / sqrt(stereoTsosGMEr.cyy());
0629         double pullGMZ_ts_stereo = (stereoTsosGMom.z() - stereoShitGMom.z()) / sqrt(stereoTsosGMEr.czz());
0630         //double pullGMX_ts_stereo = (stereoTsosGMom.x()-stereoShitGMom.x());
0631         //double pullGMY_ts_stereo = (stereoTsosGMom.y()-stereoShitGMom.y());
0632         //double pullGMZ_ts_stereo = (stereoTsosGMom.z()-stereoShitGMom.z());
0633 
0634         title.str("");
0635         title << "PullGM_X_" << subdetId << "-" << layerId << "_ts_stereo";
0636         hPullGM_X_ts_stereo[title.str()]->Fill(pullGMX_ts_stereo);
0637         title.str("");
0638         title << "PullGM_Y_" << subdetId << "-" << layerId << "_ts_stereo";
0639         hPullGM_Y_ts_stereo[title.str()]->Fill(pullGMY_ts_stereo);
0640         title.str("");
0641         title << "PullGM_Z_" << subdetId << "-" << layerId << "_ts_stereo";
0642         hPullGM_Z_ts_stereo[title.str()]->Fill(pullGMZ_ts_stereo);
0643       }
0644       lastState = updatedState;
0645       //#endif
0646     }
0647   }
0648   LogTrace("TestSmoothHits") << "end of event" << std::endl;
0649 }
0650 
0651 void TestSmoothHits::endJob() {
0652   //file->Write();
0653   TDirectory* chi2i = file->mkdir("Chi2_Increment");
0654 
0655   TDirectory* gp_ts = file->mkdir("GP_TSOS-SimHit");
0656   TDirectory* gm_ts = file->mkdir("GM_TSOS-SimHit");
0657   TDirectory* gp_tr = file->mkdir("GP_TSOS-RecHit");
0658   TDirectory* gp_rs = file->mkdir("GP_RecHit-SimHit");
0659 
0660   TDirectory* gp_tsx = gp_ts->mkdir("X");
0661   TDirectory* gp_tsy = gp_ts->mkdir("Y");
0662   TDirectory* gp_tsz = gp_ts->mkdir("Z");
0663   TDirectory* gm_tsx = gm_ts->mkdir("X");
0664   TDirectory* gm_tsy = gm_ts->mkdir("Y");
0665   TDirectory* gm_tsz = gm_ts->mkdir("Z");
0666   TDirectory* gp_trx = gp_tr->mkdir("X");
0667   TDirectory* gp_try = gp_tr->mkdir("Y");
0668   TDirectory* gp_trz = gp_tr->mkdir("Z");
0669   TDirectory* gp_rsx = gp_rs->mkdir("X");
0670   TDirectory* gp_rsy = gp_rs->mkdir("Y");
0671   TDirectory* gp_rsz = gp_rs->mkdir("Z");
0672 
0673   TDirectory* gp_tsx_mono = gp_ts->mkdir("MONOX");
0674   TDirectory* gp_tsy_mono = gp_ts->mkdir("MONOY");
0675   TDirectory* gp_tsz_mono = gp_ts->mkdir("MONOZ");
0676   TDirectory* gm_tsx_mono = gm_ts->mkdir("MONOX");
0677   TDirectory* gm_tsy_mono = gm_ts->mkdir("MONOY");
0678   TDirectory* gm_tsz_mono = gm_ts->mkdir("MONOZ");
0679   TDirectory* gp_trx_mono = gp_tr->mkdir("MONOX");
0680   TDirectory* gp_try_mono = gp_tr->mkdir("MONOY");
0681   TDirectory* gp_trz_mono = gp_tr->mkdir("MONOZ");
0682   TDirectory* gp_rsx_mono = gp_rs->mkdir("MONOX");
0683   TDirectory* gp_rsy_mono = gp_rs->mkdir("MONOY");
0684   TDirectory* gp_rsz_mono = gp_rs->mkdir("MONOZ");
0685 
0686   TDirectory* gp_tsx_stereo = gp_ts->mkdir("STEREOX");
0687   TDirectory* gp_tsy_stereo = gp_ts->mkdir("STEREOY");
0688   TDirectory* gp_tsz_stereo = gp_ts->mkdir("STEREOZ");
0689   TDirectory* gm_tsx_stereo = gm_ts->mkdir("STEREOX");
0690   TDirectory* gm_tsy_stereo = gm_ts->mkdir("STEREOY");
0691   TDirectory* gm_tsz_stereo = gm_ts->mkdir("STEREOZ");
0692   TDirectory* gp_trx_stereo = gp_tr->mkdir("STEREOX");
0693   TDirectory* gp_try_stereo = gp_tr->mkdir("STEREOY");
0694   TDirectory* gp_trz_stereo = gp_tr->mkdir("STEREOZ");
0695   TDirectory* gp_rsx_stereo = gp_rs->mkdir("STEREOX");
0696   TDirectory* gp_rsy_stereo = gp_rs->mkdir("STEREOY");
0697   TDirectory* gp_rsz_stereo = gp_rs->mkdir("STEREOZ");
0698 
0699   chi2i->cd();
0700   hTotChi2Increment->Write();
0701   hChi2_vs_Process->Write();
0702   hChi2_vs_clsize->Write();
0703   for (int i = 0; i != 6; i++)
0704     for (int j = 0; j != 9; j++) {
0705       if (i == 0 && j > 2)
0706         break;
0707       if (i == 1 && j > 1)
0708         break;
0709       if (i == 2 && j > 3)
0710         break;
0711       if (i == 3 && j > 2)
0712         break;
0713       if (i == 4 && j > 5)
0714         break;
0715       if (i == 5 && j > 8)
0716         break;
0717       chi2i->cd();
0718       title.str("");
0719       title << "Chi2Increment_" << i + 1 << "-" << j + 1;
0720       hChi2Increment[title.str()]->Write();
0721 
0722       gp_ts->cd();
0723       gp_tsx->cd();
0724       title.str("");
0725       title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_ts";
0726       hPullGP_X_ts[title.str()]->Write();
0727       gp_tsy->cd();
0728       title.str("");
0729       title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_ts";
0730       hPullGP_Y_ts[title.str()]->Write();
0731       gp_tsz->cd();
0732       title.str("");
0733       title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_ts";
0734       hPullGP_Z_ts[title.str()]->Write();
0735 
0736       gm_ts->cd();
0737       gm_tsx->cd();
0738       title.str("");
0739       title << "PullGM_X_" << i + 1 << "-" << j + 1 << "_ts";
0740       hPullGM_X_ts[title.str()]->Write();
0741       gm_tsy->cd();
0742       title.str("");
0743       title << "PullGM_Y_" << i + 1 << "-" << j + 1 << "_ts";
0744       hPullGM_Y_ts[title.str()]->Write();
0745       gm_tsz->cd();
0746       title.str("");
0747       title << "PullGM_Z_" << i + 1 << "-" << j + 1 << "_ts";
0748       hPullGM_Z_ts[title.str()]->Write();
0749 
0750       gp_tr->cd();
0751       gp_trx->cd();
0752       title.str("");
0753       title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_tr";
0754       hPullGP_X_tr[title.str()]->Write();
0755       gp_try->cd();
0756       title.str("");
0757       title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_tr";
0758       hPullGP_Y_tr[title.str()]->Write();
0759       gp_trz->cd();
0760       title.str("");
0761       title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_tr";
0762       hPullGP_Z_tr[title.str()]->Write();
0763 
0764       gp_rs->cd();
0765       gp_rsx->cd();
0766       title.str("");
0767       title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_rs";
0768       hPullGP_X_rs[title.str()]->Write();
0769       gp_rsy->cd();
0770       title.str("");
0771       title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_rs";
0772       hPullGP_Y_rs[title.str()]->Write();
0773       gp_rsz->cd();
0774       title.str("");
0775       title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_rs";
0776       hPullGP_Z_rs[title.str()]->Write();
0777 
0778       if (((i == 2 || i == 4) && (j == 0 || j == 1)) || (i == 3 || i == 5)) {
0779         //mono
0780         gp_ts->cd();
0781         gp_tsx_mono->cd();
0782         title.str("");
0783         title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_ts_mono";
0784         hPullGP_X_ts_mono[title.str()]->Write();
0785         gp_tsy_mono->cd();
0786         title.str("");
0787         title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_ts_mono";
0788         hPullGP_Y_ts_mono[title.str()]->Write();
0789         gp_tsz_mono->cd();
0790         title.str("");
0791         title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_ts_mono";
0792         hPullGP_Z_ts_mono[title.str()]->Write();
0793 
0794         gm_ts->cd();
0795         gm_tsx_mono->cd();
0796         title.str("");
0797         title << "PullGM_X_" << i + 1 << "-" << j + 1 << "_ts_mono";
0798         hPullGM_X_ts_mono[title.str()]->Write();
0799         gm_tsy_mono->cd();
0800         title.str("");
0801         title << "PullGM_Y_" << i + 1 << "-" << j + 1 << "_ts_mono";
0802         hPullGM_Y_ts_mono[title.str()]->Write();
0803         gm_tsz_mono->cd();
0804         title.str("");
0805         title << "PullGM_Z_" << i + 1 << "-" << j + 1 << "_ts_mono";
0806         hPullGM_Z_ts_mono[title.str()]->Write();
0807 
0808         gp_tr->cd();
0809         gp_trx_mono->cd();
0810         title.str("");
0811         title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_tr_mono";
0812         hPullGP_X_tr_mono[title.str()]->Write();
0813         gp_try_mono->cd();
0814         title.str("");
0815         title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_tr_mono";
0816         hPullGP_Y_tr_mono[title.str()]->Write();
0817         gp_trz_mono->cd();
0818         title.str("");
0819         title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_tr_mono";
0820         hPullGP_Z_tr_mono[title.str()]->Write();
0821 
0822         gp_rs->cd();
0823         gp_rsx_mono->cd();
0824         title.str("");
0825         title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_rs_mono";
0826         hPullGP_X_rs_mono[title.str()]->Write();
0827         gp_rsy_mono->cd();
0828         title.str("");
0829         title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_rs_mono";
0830         hPullGP_Y_rs_mono[title.str()]->Write();
0831         gp_rsz_mono->cd();
0832         title.str("");
0833         title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_rs_mono";
0834         hPullGP_Z_rs_mono[title.str()]->Write();
0835 
0836         //stereo
0837         gp_ts->cd();
0838         gp_tsx_stereo->cd();
0839         title.str("");
0840         title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0841         hPullGP_X_ts_stereo[title.str()]->Write();
0842         gp_tsy_stereo->cd();
0843         title.str("");
0844         title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0845         hPullGP_Y_ts_stereo[title.str()]->Write();
0846         gp_tsz_stereo->cd();
0847         title.str("");
0848         title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0849         hPullGP_Z_ts_stereo[title.str()]->Write();
0850 
0851         gm_ts->cd();
0852         gm_tsx_stereo->cd();
0853         title.str("");
0854         title << "PullGM_X_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0855         hPullGM_X_ts_stereo[title.str()]->Write();
0856         gm_tsy_stereo->cd();
0857         title.str("");
0858         title << "PullGM_Y_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0859         hPullGM_Y_ts_stereo[title.str()]->Write();
0860         gm_tsz_stereo->cd();
0861         title.str("");
0862         title << "PullGM_Z_" << i + 1 << "-" << j + 1 << "_ts_stereo";
0863         hPullGM_Z_ts_stereo[title.str()]->Write();
0864 
0865         gp_tr->cd();
0866         gp_trx_stereo->cd();
0867         title.str("");
0868         title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_tr_stereo";
0869         hPullGP_X_tr_stereo[title.str()]->Write();
0870         gp_try_stereo->cd();
0871         title.str("");
0872         title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_tr_stereo";
0873         hPullGP_Y_tr_stereo[title.str()]->Write();
0874         gp_trz_stereo->cd();
0875         title.str("");
0876         title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_tr_stereo";
0877         hPullGP_Z_tr_stereo[title.str()]->Write();
0878 
0879         gp_rs->cd();
0880         gp_rsx_stereo->cd();
0881         title.str("");
0882         title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_rs_stereo";
0883         hPullGP_X_rs_stereo[title.str()]->Write();
0884         gp_rsy_stereo->cd();
0885         title.str("");
0886         title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_rs_stereo";
0887         hPullGP_Y_rs_stereo[title.str()]->Write();
0888         gp_rsz_stereo->cd();
0889         title.str("");
0890         title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_rs_stereo";
0891         hPullGP_Z_rs_stereo[title.str()]->Write();
0892       }
0893     }
0894 
0895   file->Close();
0896 }
0897 
0898 //needed by to do the residual for matched hits
0899 //taken from SiStripTrackingRecHitsValid.cc
0900 std::pair<LocalPoint, LocalVector> TestSmoothHits::projectHit(const PSimHit& hit,
0901                                                               const StripGeomDetUnit* stripDet,
0902                                                               const BoundPlane& plane) {
0903   const StripTopology& topol = stripDet->specificTopology();
0904   GlobalPoint globalpos = stripDet->surface().toGlobal(hit.localPosition());
0905   LocalPoint localHit = plane.toLocal(globalpos);
0906   //track direction
0907   LocalVector locdir = hit.localDirection();
0908   //rotate track in new frame
0909 
0910   GlobalVector globaldir = stripDet->surface().toGlobal(locdir);
0911   LocalVector dir = plane.toLocal(globaldir);
0912   float scale = -localHit.z() / dir.z();
0913 
0914   LocalPoint projectedPos = localHit + scale * dir;
0915 
0916   float selfAngle = topol.stripAngle(topol.strip(hit.localPosition()));
0917 
0918   LocalVector stripDir(sin(selfAngle), cos(selfAngle), 0);  // vector along strip in hit frame
0919 
0920   LocalVector localStripDir(plane.toLocal(stripDet->surface().toGlobal(stripDir)));
0921 
0922   return std::pair<LocalPoint, LocalVector>(projectedPos, localStripDir);
0923 }
0924 
0925 #include "FWCore/Framework/interface/ModuleFactory.h"
0926 #include "FWCore/Framework/interface/MakerMacros.h"
0927 
0928 DEFINE_FWK_MODULE(TestSmoothHits);