Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:56 

0001 #include "CkfDebugger.h"
0002 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0003 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0004 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0005 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0006 #include "Geometry/CommonTopologies/interface/Topology.h"
0007 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementVector.h"
0008 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0009 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0010 #include "TrackingTools/DetLayers/interface/NavigationSchool.h"
0011 #include "RecoTracker/Record/interface/NavigationSchoolRecord.h"
0012 #include "TSOSFromSimHitFactory.h"
0013 #include "TrackingTools/MeasurementDet/interface/MeasurementDet.h"
0014 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0015 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0016 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0017 
0018 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0019 
0020 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0021 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0022 #include "TrackingTools/DetLayers/interface/TkLayerLess.h"
0023 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0024 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0025 // #include "RecoTracker/TkDetLayers/interface/PixelBarrelLayer.h"
0026 
0027 #include "TrackingTools/TransientTrackingRecHit/interface/TrackingRecHitProjector.h"
0028 #include "RecoTracker/TransientTrackingRecHit/interface/ProjectedRecHit2D.h"
0029 
0030 #include <iostream>
0031 #include <iomanip>
0032 #include <sstream>
0033 
0034 using namespace std;
0035 
0036 inline DetId gluedId(const DetId& du) {
0037   unsigned int mask = ~3;  // mask the last two bits
0038   return DetId(du.rawId() & mask);
0039 }
0040 
0041 CkfDebugger::CkfDebugger(edm::ConsumesCollector iC)
0042     : theTrackerToken(iC.esConsumes()),
0043       theFieldToken(iC.esConsumes()),
0044       theTopoHandToken(iC.esConsumes()),
0045       theNavToken(iC.esConsumes(edm::ESInputTag("", "SimpleNavigationSchool"))),
0046       trackerHitAssociatorConfig_(std::move(iC)),
0047       totSeeds(0) {}
0048 
0049 void CkfDebugger::setConditions(edm::EventSetup const& es) {
0050   file = new TFile("out.root", "recreate");
0051   hchi2seedAll = new TH1F("hchi2seedAll", "hchi2seedAll", 2000, 0, 200);
0052   hchi2seedProb = new TH1F("hchi2seedProb", "hchi2seedProb", 2000, 0, 200);
0053 
0054   theTrackerGeom = &es.getData(theTrackerToken);
0055 
0056   theMagField = &es.getData(theFieldToken);
0057 
0058   //Retrieve tracker topology from geometry
0059   theTopo = &es.getData(theTopoHandToken);
0060 
0061   theNavSchool = &es.getData(theNavToken);
0062 
0063   for (int i = 0; i != 17; i++) {
0064     dump.push_back(0);
0065   }
0066 
0067   std::stringstream title;
0068   for (int i = 0; i != 6; i++)
0069     for (int j = 0; j != 9; j++) {
0070       if (i == 0 && j > 2)
0071         break;
0072       if (i == 1 && j > 1)
0073         break;
0074       if (i == 2 && j > 3)
0075         break;
0076       if (i == 3 && j > 2)
0077         break;
0078       if (i == 4 && j > 5)
0079         break;
0080       if (i == 5 && j > 8)
0081         break;
0082       dump2[pair<int, int>(i, j)] = 0;
0083       dump3[pair<int, int>(i, j)] = 0;
0084       dump4[pair<int, int>(i, j)] = 0;
0085       dump5[pair<int, int>(i, j)] = 0;
0086       dump6[pair<int, int>(i, j)] = 0;
0087       title.str("");
0088       title << "pullX_" << i + 1 << "-" << j + 1 << "_sh-rh";
0089       hPullX_shrh[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0090       title.str("");
0091       title << "pullY_" << i + 1 << "-" << j + 1 << "_sh-rh";
0092       hPullY_shrh[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0093       title.str("");
0094       title << "pullX_" << i + 1 << "-" << j + 1 << "_sh-st";
0095       hPullX_shst[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0096       title.str("");
0097       title << "pullY_" << i + 1 << "-" << j + 1 << "_sh-st";
0098       hPullY_shst[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0099       title.str("");
0100       title << "pullX_" << i + 1 << "-" << j + 1 << "_st-rh";
0101       hPullX_strh[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0102       title.str("");
0103       title << "pullY_" << i + 1 << "-" << j + 1 << "_st-rh";
0104       hPullY_strh[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0105       title.str("");
0106       title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_sh-st";
0107       hPullGP_X_shst[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0108       title.str("");
0109       title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_sh-st";
0110       hPullGP_Y_shst[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0111       title.str("");
0112       title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_sh-st";
0113       hPullGP_Z_shst[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0114       if (((i == 2 || i == 4) && (j == 0 || j == 1)) || (i == 3 || i == 5)) {
0115         title.str("");
0116         title << "pullM_" << i + 1 << "-" << j + 1 << "_sh-rh";
0117         hPullM_shrh[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0118         title.str("");
0119         title << "pullS_" << i + 1 << "-" << j + 1 << "_sh-rh";
0120         hPullS_shrh[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0121         title.str("");
0122         title << "pullM_" << i + 1 << "-" << j + 1 << "_sh-st";
0123         hPullM_shst[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0124         title.str("");
0125         title << "pullS_" << i + 1 << "-" << j + 1 << "_sh-st";
0126         hPullS_shst[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0127         title.str("");
0128         title << "pullM_" << i + 1 << "-" << j + 1 << "_st-rh";
0129         hPullM_strh[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0130         title.str("");
0131         title << "pullS_" << i + 1 << "-" << j + 1 << "_st-rh";
0132         hPullS_strh[title.str()] = new TH1F(title.str().c_str(), title.str().c_str(), 1000, -50, 50);
0133       }
0134     }
0135 
0136   hPullGPXvsGPX_shst = new TH2F("PullGPXvsGPX_shst", "PullGPXvsGPX_shst", 1000, -50, 50, 100, -50, 50);
0137   hPullGPXvsGPY_shst = new TH2F("PullGPXvsGPY_shst", "PullGPXvsGPY_shst", 1000, -50, 50, 100, -50, 50);
0138   hPullGPXvsGPZ_shst = new TH2F("PullGPXvsGPZ_shst", "PullGPXvsGPZ_shst", 1000, -50, 50, 200, -100, 100);
0139   hPullGPXvsGPr_shst = new TH2F("PullGPXvsGPr_shst", "PullGPXvsGPr_shst", 1000, -50, 50, 300, -150, 150);
0140   hPullGPXvsGPeta_shst = new TH2F("PullGPXvsGPeta_shst", "PullGPXvsGPeta_shst", 1000, -50, 50, 50, -2.5, 2.5);
0141   hPullGPXvsGPphi_shst = new TH2F("PullGPXvsGPphi_shst", "PullGPXvsGPphi_shst", 1000, -50, 50, 63, 0, 6.3);
0142 
0143   seedWithDelta = 0;
0144   problems = 0;
0145   no_sim_hit = 0;
0146   no_layer = 0;
0147   layer_not_found = 0;
0148   det_not_found = 0;
0149   chi2gt30 = 0;
0150   chi2gt30delta = 0;
0151   chi2gt30deltaSeed = 0;
0152   chi2ls30 = 0;
0153   simple_hit_not_found = 0;
0154   no_component = 0;
0155   only_one_component = 0;
0156   matched_not_found = 0;
0157   matched_not_associated = 0;
0158   partner_det_not_fuond = 0;
0159   glued_det_not_fuond = 0;
0160   propagation = 0;
0161   other = 0;
0162   totchi2gt30 = 0;
0163 }
0164 
0165 void CkfDebugger::printSimHits(const edm::Event& iEvent) {
0166   edm::LogVerbatim("CkfDebugger") << "\nEVENT #" << iEvent.id();
0167 
0168   hitAssociator = new TrackerHitAssociator(
0169       iEvent, trackerHitAssociatorConfig_);  //delete deleteHitAssociator() in TrackCandMaker.cc
0170 
0171   std::map<unsigned int, std::vector<PSimHit> >& theHitsMap = hitAssociator->SimHitMap;
0172   idHitsMap.clear();
0173 
0174   for (std::map<unsigned int, std::vector<PSimHit> >::iterator it = theHitsMap.begin(); it != theHitsMap.end(); it++) {
0175     for (std::vector<PSimHit>::iterator isim = it->second.begin(); isim != it->second.end(); ++isim) {
0176       idHitsMap[isim->trackId()].push_back(&*isim);
0177     }
0178   }
0179 
0180   for (std::map<unsigned int, std::vector<PSimHit*> >::iterator it = idHitsMap.begin(); it != idHitsMap.end(); it++) {
0181     sort(it->second.begin(), it->second.end(), [](auto* a, auto* b) { return a->timeOfFlight() < b->timeOfFlight(); });
0182     for (std::vector<PSimHit*>::iterator isim = it->second.begin(); isim != it->second.end(); ++isim) {
0183       const GeomDetUnit* detUnit = theTrackerGeom->idToDetUnit(DetId((*isim)->detUnitId()));
0184       dumpSimHit(SimHit((*isim), detUnit));
0185     }
0186   }
0187 }
0188 
0189 void CkfDebugger::dumpSimHit(const SimHit& hit) const {
0190   GlobalPoint pos = hit.globalPosition();
0191   edm::LogVerbatim("CkfDebugger") << "SimHit pos" << pos << " r=" << pos.perp() << " phi=" << pos.phi()
0192                                   << " trackId=" << hit.trackId() << " particleType=" << hit.particleType()
0193                                   << " pabs=" << hit.pabs() << " processType=" << hit.processType();
0194 }
0195 
0196 bool CkfDebugger::analyseCompatibleMeasurements(const Trajectory& traj,
0197                                                 const std::vector<TrajectoryMeasurement>& meas,
0198                                                 const MeasurementTrackerEvent* aMeasurementTracker,
0199                                                 const Propagator* propagator,
0200                                                 const Chi2MeasurementEstimatorBase* estimator,
0201                                                 const TransientTrackingRecHitBuilder* aTTRHBuilder) {
0202   LogTrace("CkfDebugger") << "\nnow in analyseCompatibleMeasurements";
0203   LogTrace("CkfDebugger") << "number of input hits:" << meas.size();
0204   for (std::vector<TrajectoryMeasurement>::const_iterator tmpIt = meas.begin(); tmpIt != meas.end(); tmpIt++) {
0205     if (tmpIt->recHit()->isValid())
0206       LogTrace("CkfDebugger") << "valid hit at position:" << tmpIt->recHit()->globalPosition();
0207   }
0208   theForwardPropagator = propagator;
0209   theChi2 = estimator;
0210   theMeasurementTracker = aMeasurementTracker;
0211   theGeomSearchTracker = theMeasurementTracker->geometricSearchTracker();
0212   theTTRHBuilder = aTTRHBuilder;
0213   unsigned int trajId = 0;
0214   if (!correctTrajectory(traj, trajId)) {
0215     LogTrace("CkfDebugger") << "trajectory not correct";
0216     return true;
0217   }  // only correct trajectories analysed
0218   LogTrace("CkfDebugger") << "correct trajectory";
0219 
0220   if (traj.measurements().size() == 2) {
0221     if (testSeed(traj.firstMeasurement().recHit(),
0222                  traj.lastMeasurement().recHit(),
0223                  traj.lastMeasurement().updatedState()) == -1) {
0224       LogTrace("CkfDebugger") << "Seed has delta";
0225       seedWithDelta++;
0226       return false;  //true;//false?
0227     }
0228   }
0229 
0230   //const PSimHit* correctHit = nextCorrectHit(traj, trajId);
0231   //if ( correctHit == 0) return true; // no more simhits on this track
0232   std::vector<const PSimHit*> correctHits = nextCorrectHits(traj, trajId);
0233   if (correctHits.empty())
0234     return true;  // no more simhits on this track
0235 
0236   for (std::vector<const PSimHit*>::iterator corHit = correctHits.begin(); corHit != correctHits.end(); corHit++) {
0237     for (std::vector<TM>::const_iterator i = meas.begin(); i != meas.end(); i++) {
0238       if (correctMeas(*i, *corHit)) {
0239         LogTrace("CkfDebugger") << "Correct hit found at position " << i - meas.begin();
0240         return true;
0241       }
0242     }
0243   }
0244 
0245   //debug why the first hit in correctHits is not found
0246   //FIXME should loop over all hits
0247   const PSimHit* correctHit = *(correctHits.begin());
0248 
0249   // correct hit not found
0250   edm::LogVerbatim("CkfDebugger") << std::endl
0251                                   << "CkfDebugger: problem found: correct hit not found by findCompatibleMeasurements";
0252   edm::LogVerbatim("CkfDebugger") << "The correct hit position is " << position(correctHit) << " lp "
0253                                   << correctHit->localPosition();
0254   edm::LogVerbatim("CkfDebugger") << "The size of the meas vector is " << meas.size();
0255   dump[0]++;
0256   problems++;
0257 
0258   for (std::vector<TM>::const_iterator i = meas.begin(); i != meas.end(); i++) {
0259     edm::LogVerbatim("CkfDebugger") << "Is the hit valid? " << i->recHit()->isValid();
0260     if (i->recHit()->isValid()) {
0261       edm::LogVerbatim("CkfDebugger") << "RecHit at " << i->recHit()->globalPosition() << " layer "
0262                                       << ((i->recHit()->det()->geographicalId().rawId() >> 16) & 0xF) << " subdet "
0263                                       << i->recHit()->det()->geographicalId().subdetId() << " Chi2 " << i->estimate();
0264     } else if (i->recHit()->det() == nullptr) {
0265       edm::LogVerbatim("CkfDebugger") << "Invalid RecHit returned with zero Det pointer";
0266     } else if (i->recHit()->det() == det(correctHit)) {
0267       edm::LogVerbatim("CkfDebugger") << "Invalid hit returned in correct Det";
0268     } else {
0269       edm::LogVerbatim("CkfDebugger") << "Invalid hit returned in Det at gpos " << i->recHit()->det()->position()
0270                                       << " correct Det is at " << det(correctHit)->position();
0271     }
0272   }
0273 
0274   //Look if the correct RecHit exists
0275   std::pair<CTTRHp, double> correctRecHit = analyseRecHitExistance(*correctHit, traj.lastMeasurement().updatedState());
0276   if (correctRecHit.first == nullptr) {
0277     //the hit does not exist or is uncorrectly matched
0278     if (fabs(correctRecHit.second - 0) < 0.01) {
0279       dump[1]++;
0280     }  //other
0281     if (fabs(correctRecHit.second + 1) < 0.01) {
0282       dump[8]++;
0283     }  //propagation
0284     if (fabs(correctRecHit.second + 2) < 0.01) {
0285       dump[9]++;
0286     }  //No component is found
0287     if (fabs(correctRecHit.second + 3) < 0.01) {
0288       dump[10]++;
0289     }  //Partner measurementDet not found
0290     if (fabs(correctRecHit.second + 4) < 0.01) {
0291       dump[11]++;
0292     }  //glued MeasurementDet not found
0293     if (fabs(correctRecHit.second + 5) < 0.01) {
0294       dump[12]++;
0295     }  //matched not found
0296     if (fabs(correctRecHit.second + 6) < 0.01) {
0297       dump[13]++;
0298     }  //Matched not associated
0299     if (fabs(correctRecHit.second + 7) < 0.01) {
0300       dump[14]++;
0301     }  //Only one component is found
0302     if (fabs(correctRecHit.second + 8) < 0.01) {
0303       dump[15]++;
0304     }  //not found (is not a glued det)
0305   } else {
0306     //the hit exists: why wasn't it found?
0307     int result = analyseRecHitNotFound(traj, correctRecHit.first);
0308     if (result == 5) {
0309       if (correctRecHit.second > 30) {
0310         edm::LogVerbatim("CkfDebugger") << "Outling RecHit at pos=" << correctRecHit.first->globalPosition()
0311                                         << " from SimHit at pos=" << position(correctHit)
0312                                         << " det=" << correctHit->detUnitId()
0313                                         << " process=" << correctHit->processType();
0314         if (hasDelta(correctHit)) {
0315           edm::LogVerbatim("CkfDebugger") << "there are deltas on this det";
0316           chi2gt30delta++;
0317           dump5[pair<int, int>((correctRecHit.first->det()->geographicalId().subdetId() - 1),
0318                                (layer(correctRecHit.first->det())) - 1)]++;
0319         } else {
0320           edm::LogVerbatim("CkfDebugger") << "no deltas on this det";
0321           dump[5]++;
0322           chi2gt30++;
0323           dump3[pair<int, int>((correctRecHit.first->det()->geographicalId().subdetId() - 1),
0324                                (layer(correctRecHit.first->det())) - 1)]++;
0325           CTTRHp h1 = traj.measurements()[0].recHit();
0326           CTTRHp h2 = traj.measurements()[1].recHit();
0327           TSOS t = traj.measurements()[1].updatedState();
0328           double chi2 = testSeed(h1, h2, t);
0329           if (chi2 == -1) {
0330             edm::LogVerbatim("CkfDebugger") << "there were deltas in the seed";
0331             chi2gt30deltaSeed++;
0332           } else {
0333             hchi2seedProb->Fill(chi2);
0334             edm::LogVerbatim("CkfDebugger") << "no deltas in the seed. What is wrong?";
0335 
0336             TSOS detState = theForwardPropagator->propagate(traj.lastMeasurement().updatedState(),
0337                                                             correctRecHit.first->det()->surface());
0338             TSOS simDetState =
0339                 theForwardPropagator->propagate(traj.lastMeasurement().updatedState(), det(correctHit)->surface());
0340 
0341             if (true /*detState.globalMomentum().y()>0*/) {
0342               int subdetId = correctRecHit.first->det()->geographicalId().subdetId();
0343               int layerId = layer(correctRecHit.first->det());
0344 
0345               LogTrace("CkfDebugger") << "position(correctHit)=" << position(correctHit);
0346               LogTrace("CkfDebugger") << "correctRecHit.first->globalPosition()="
0347                                       << correctRecHit.first->globalPosition();
0348               LogTrace("CkfDebugger") << "detState.globalPosition()=" << detState.globalPosition();
0349               LogTrace("CkfDebugger") << "simDetState.globalPosition()=" << simDetState.globalPosition();
0350 
0351               LogTrace("CkfDebugger") << "correctHit->localPosition()=" << correctHit->localPosition();
0352               LogTrace("CkfDebugger") << "correctRecHit.first->localPosition()="
0353                                       << correctRecHit.first->localPosition();
0354               LogTrace("CkfDebugger") << "correctRecHit.first->localPositionError()="
0355                                       << correctRecHit.first->localPositionError();
0356               LogTrace("CkfDebugger") << "detState.localPosition()=" << detState.localPosition();
0357               LogTrace("CkfDebugger") << "detState.localError().positionError()="
0358                                       << detState.localError().positionError();
0359               LogTrace("CkfDebugger") << "simDetState.localPosition()=" << simDetState.localPosition();
0360               LogTrace("CkfDebugger") << "simDetState.localError().positionError()="
0361                                       << simDetState.localError().positionError();
0362               double pullx_shrh = (correctHit->localPosition().x() - correctRecHit.first->localPosition().x()) /
0363                                   sqrt(correctRecHit.first->localPositionError().xx());
0364               double pully_shrh = 0;
0365               if (correctRecHit.first->localPositionError().yy() != 0)
0366                 pully_shrh = (correctHit->localPosition().y() - correctRecHit.first->localPosition().y()) /
0367                              sqrt(correctRecHit.first->localPositionError().yy());
0368               double pullx_shst = (correctHit->localPosition().x() - simDetState.localPosition().x()) /
0369                                   sqrt(simDetState.localError().positionError().xx());
0370               double pully_shst = (correctHit->localPosition().y() - simDetState.localPosition().y()) /
0371                                   sqrt(simDetState.localError().positionError().yy());
0372 
0373               LogTrace("CkfDebugger") << "pullx(sh-rh)=" << pullx_shrh;
0374               LogTrace("CkfDebugger") << "pully(sh-rh)=" << pully_shrh;
0375               LogTrace("CkfDebugger") << "pullx(sh-st)=" << pullx_shst;
0376               LogTrace("CkfDebugger") << "pully(sh-st)=" << pully_shst;
0377 
0378               LogTrace("CkfDebugger") << "pullx(st-rh)="
0379                                       << (detState.localPosition().x() - correctRecHit.first->localPosition().x()) /
0380                                              sqrt(correctRecHit.first->localPositionError().xx() +
0381                                                   detState.localError().positionError().xx());
0382 
0383               std::pair<double, double> pulls = computePulls(correctRecHit.first, detState);
0384               if (subdetId > 0 && subdetId < 7 && layerId > 0 && layerId < 10) {
0385                 stringstream title;
0386                 title.str("");
0387                 title << "pullX_" << subdetId << "-" << layerId << "_sh-rh";
0388                 hPullX_shrh[title.str()]->Fill(pullx_shrh);
0389                 title.str("");
0390                 title << "pullY_" << subdetId << "-" << layerId << "_sh-rh";
0391                 hPullY_shrh[title.str()]->Fill(pully_shrh);
0392                 title.str("");
0393                 title << "pullX_" << subdetId << "-" << layerId << "_sh-st";
0394                 hPullX_shst[title.str()]->Fill(pullx_shst);
0395                 title.str("");
0396                 title << "pullY_" << subdetId << "-" << layerId << "_sh-st";
0397                 hPullY_shst[title.str()]->Fill(pully_shst);
0398                 title.str("");
0399                 title << "pullX_" << subdetId << "-" << layerId << "_st-rh";
0400                 hPullX_strh[title.str()]->Fill(pulls.first);
0401                 title.str("");
0402                 title << "pullY_" << subdetId << "-" << layerId << "_st-rh";
0403                 hPullY_strh[title.str()]->Fill(pulls.second);
0404 
0405                 GlobalPoint shGPos = position(correctHit);
0406                 GlobalPoint stGPos = simDetState.globalPosition();
0407                 GlobalError stGPosErr = simDetState.cartesianError().position();
0408                 double pullGPx = (shGPos.x() - stGPos.x()) / sqrt(stGPosErr.cxx());
0409                 title.str("");
0410                 title << "PullGP_X_" << subdetId << "-" << layerId << "_sh-st";
0411                 hPullGP_X_shst[title.str()]->Fill(pullGPx);
0412                 title.str("");
0413                 title << "PullGP_Y_" << subdetId << "-" << layerId << "_sh-st";
0414                 hPullGP_Y_shst[title.str()]->Fill((shGPos.y() - stGPos.y()) / sqrt(stGPosErr.cyy()));
0415                 title.str("");
0416                 title << "PullGP_Z_" << subdetId << "-" << layerId << "_sh-st";
0417                 hPullGP_Z_shst[title.str()]->Fill((shGPos.z() - stGPos.z()) / sqrt(stGPosErr.czz()));
0418 
0419                 if (subdetId == 3 && layerId == 1) {
0420                   hPullGPXvsGPX_shst->Fill(pullGPx, shGPos.x());
0421                   hPullGPXvsGPY_shst->Fill(pullGPx, shGPos.y());
0422                   hPullGPXvsGPZ_shst->Fill(pullGPx, shGPos.z());
0423                   hPullGPXvsGPr_shst->Fill(pullGPx, shGPos.mag());
0424                   hPullGPXvsGPeta_shst->Fill(pullGPx, shGPos.eta());
0425                   hPullGPXvsGPphi_shst->Fill(pullGPx, shGPos.phi());
0426                 }
0427                 if (dynamic_cast<const SiStripMatchedRecHit2D*>(correctRecHit.first->hit())) {
0428                   LogTrace("CkfDebugger") << "MONO HIT";
0429                   auto m = dynamic_cast<const SiStripMatchedRecHit2D*>(correctRecHit.first->hit())->monoHit();
0430                   CTTRHp tMonoHit = theTTRHBuilder->build(&m);
0431                   const PSimHit sMonoHit = *(hitAssociator->associateHit(*tMonoHit->hit()).begin());
0432                   TSOS monoState = theForwardPropagator->propagate(traj.lastMeasurement().updatedState(),
0433                                                                    tMonoHit->det()->surface());
0434                   double pullM_shrh = (sMonoHit.localPosition().x() - tMonoHit->localPosition().x()) /
0435                                       sqrt(tMonoHit->localPositionError().xx());
0436                   double pullM_shst = (sMonoHit.localPosition().x() - monoState.localPosition().x()) /
0437                                       sqrt(monoState.localError().positionError().xx());
0438                   std::pair<double, double> pullsMono = computePulls(tMonoHit, monoState);
0439                   title.str("");
0440                   title << "pullM_" << subdetId << "-" << layerId << "_sh-rh";
0441                   hPullM_shrh[title.str()]->Fill(pullM_shrh);
0442                   title.str("");
0443                   title << "pullM_" << subdetId << "-" << layerId << "_sh-st";
0444                   hPullM_shst[title.str()]->Fill(pullM_shst);
0445                   title.str("");
0446                   title << "pullM_" << subdetId << "-" << layerId << "_st-rh";
0447                   hPullM_strh[title.str()]->Fill(pullsMono.first);
0448 
0449                   LogTrace("CkfDebugger") << "STEREO HIT";
0450                   auto s = dynamic_cast<const SiStripMatchedRecHit2D*>(correctRecHit.first->hit())->stereoHit();
0451                   CTTRHp tStereoHit = theTTRHBuilder->build(&s);
0452                   const PSimHit sStereoHit = *(hitAssociator->associateHit(*tStereoHit->hit()).begin());
0453                   TSOS stereoState = theForwardPropagator->propagate(traj.lastMeasurement().updatedState(),
0454                                                                      tStereoHit->det()->surface());
0455                   double pullS_shrh = (sStereoHit.localPosition().x() - tStereoHit->localPosition().x()) /
0456                                       sqrt(tStereoHit->localPositionError().xx());
0457                   double pullS_shst = (sStereoHit.localPosition().x() - stereoState.localPosition().x()) /
0458                                       sqrt(stereoState.localError().positionError().xx());
0459                   std::pair<double, double> pullsStereo = computePulls(tStereoHit, stereoState);
0460                   title.str("");
0461                   title << "pullS_" << subdetId << "-" << layerId << "_sh-rh";
0462                   hPullS_shrh[title.str()]->Fill(pullS_shrh);
0463                   title.str("");
0464                   title << "pullS_" << subdetId << "-" << layerId << "_sh-st";
0465                   hPullS_shst[title.str()]->Fill(pullS_shst);
0466                   title.str("");
0467                   title << "pullS_" << subdetId << "-" << layerId << "_st-rh";
0468                   hPullS_strh[title.str()]->Fill(pullsStereo.first);
0469                 }
0470               } else
0471                 edm::LogVerbatim("CkfDebugger")
0472                     << "unexpected result: wrong det or layer id " << subdetId << " " << layerId << " "
0473                     << correctRecHit.first->det()->geographicalId().rawId();
0474             }
0475           }
0476         }
0477       } else {
0478         edm::LogVerbatim("CkfDebugger") << "unexpected result " << correctRecHit.second;
0479         dump[6]++;
0480         chi2ls30++;
0481       }
0482     } else
0483       dump[result]++;
0484     if (result == 3) {
0485       dump2[pair<int, int>((correctRecHit.first->det()->geographicalId().subdetId() - 1),
0486                            (layer(correctRecHit.first->det())) - 1)]++;
0487     }
0488     if (result == 4) {
0489       dump4[pair<int, int>((correctRecHit.first->det()->geographicalId().subdetId() - 1),
0490                            (layer(correctRecHit.first->det())) - 1)]++;
0491     }
0492     if (correctRecHit.second > 30) {
0493       dump[7]++;
0494       totchi2gt30++;
0495     }
0496   }
0497   return false;
0498 }
0499 
0500 bool CkfDebugger::correctTrajectory(const Trajectory& traj, unsigned int& trajId) const {
0501   LogTrace("CkfDebugger") << "now in correctTrajectory";
0502   Trajectory::RecHitContainer hits = traj.recHits();
0503 
0504   std::vector<SimHitIdpr> currentTrackId = hitAssociator->associateHitId(*hits.front()->hit());
0505   if (currentTrackId.empty())
0506     return false;
0507 
0508   for (Trajectory::RecHitContainer::const_iterator rh = hits.begin(); rh != hits.end(); ++rh) {
0509     //if invalid hit exit
0510     if (!(*rh)->hit()->isValid()) {
0511       //LogTrace("CkfDebugger") << "invalid hit" ;
0512       return false;
0513     }
0514 
0515     //if hits from deltas exit
0516     bool nogoodhit = true;
0517     std::vector<PSimHit> assSimHits = hitAssociator->associateHit(*(*rh)->hit());
0518     for (std::vector<PSimHit>::iterator shit = assSimHits.begin(); shit != assSimHits.end(); shit++) {
0519       if (goodSimHit(*shit))
0520         nogoodhit = false;
0521     }
0522     if (nogoodhit)
0523       return false;
0524 
0525     //all hits must be associated to the same sim track
0526     bool test = true;
0527     std::vector<SimHitIdpr> nextTrackId = hitAssociator->associateHitId(*(*rh)->hit());
0528     for (std::vector<SimHitIdpr>::iterator i = currentTrackId.begin(); i != currentTrackId.end(); i++) {
0529       for (std::vector<SimHitIdpr>::iterator j = nextTrackId.begin(); j != nextTrackId.end(); j++) {
0530         if (i->first == j->first)
0531           test = false;
0532         //LogTrace("CkfDebugger") << "valid " << *i << " " << *j ;
0533         trajId = j->first;
0534       }
0535     }
0536     if (test) { /*LogTrace("CkfDebugger") << "returning false" ;*/
0537       return false;
0538     }
0539     //     std::vector<PSimHit*> simTrackHits = idHitsMap[trajId];
0540     //     if (!goodSimHit(simTrackHits.))
0541   }
0542   //LogTrace("CkfDebugger") << "returning true" ;
0543   return true;
0544 }
0545 
0546 int CkfDebugger::assocTrackId(CTTRHp rechit) const {
0547   LogTrace("CkfDebugger") << "now in assocTrackId";
0548 
0549   if (!rechit->hit()->isValid()) {
0550     return -1;
0551   }
0552 
0553   std::vector<SimHitIdpr> ids = hitAssociator->associateHitId(*rechit->hit());
0554   if (!ids.empty()) {
0555     return ids[0].first;  //FIXME if size>1!!
0556   } else {
0557     return -1;
0558   }
0559 }
0560 
0561 vector<const PSimHit*> CkfDebugger::nextCorrectHits(const Trajectory& traj, unsigned int& trajId) {
0562   std::vector<const PSimHit*> result;
0563   // find the component of the RecHit at largest distance from origin (FIXME: should depend on propagation direction)
0564   LogTrace("CkfDebugger") << "now in nextCorrectHits";
0565   TransientTrackingRecHit::ConstRecHitPointer lastRecHit = traj.lastMeasurement().recHit();
0566   TransientTrackingRecHit::RecHitContainer comp = lastRecHit->transientHits();
0567   if (!comp.empty()) {
0568     float maxR = 0;
0569     for (TransientTrackingRecHit::RecHitContainer::const_iterator ch = comp.begin(); ch != comp.end(); ++ch) {
0570       if ((*ch)->globalPosition().mag() > maxR)
0571         lastRecHit = *ch;
0572       maxR = (*ch)->globalPosition().mag();
0573     }
0574   }
0575   edm::LogVerbatim("CkfDebugger") << "CkfDebugger: lastRecHit is at gpos " << lastRecHit->globalPosition() << " layer "
0576                                   << layer((lastRecHit->det())) << " subdet "
0577                                   << lastRecHit->det()->geographicalId().subdetId();
0578 
0579   //find the simHits associated to the recHit
0580   const std::vector<PSimHit>& pSimHitVec = hitAssociator->associateHit(*lastRecHit->hit());
0581   for (std::vector<PSimHit>::const_iterator shit = pSimHitVec.begin(); shit != pSimHitVec.end(); shit++) {
0582     const GeomDetUnit* detUnit = theTrackerGeom->idToDetUnit(DetId(shit->detUnitId()));
0583     LogTrace("CkfDebugger") << "from hitAssociator SimHits are at GP=" << detUnit->toGlobal(shit->localPosition())
0584                             << " traId=" << shit->trackId() << " particleType " << shit->particleType()
0585                             << " pabs=" << shit->pabs() << " detUnitId=" << shit->detUnitId() << " layer "
0586                             << layer((det(&*shit))) << " subdet " << det(&*shit)->geographicalId().subdetId();
0587   }
0588 
0589   //choose the simHit from the same track that has the highest tof
0590   const PSimHit* lastPSH = nullptr;
0591   if (!pSimHitVec.empty()) {
0592     float maxTOF = 0;
0593     for (std::vector<PSimHit>::const_iterator ch = pSimHitVec.begin(); ch != pSimHitVec.end(); ++ch) {
0594       if ((ch->trackId() == trajId) && (ch->timeOfFlight() > maxTOF) && (goodSimHit(*ch))) {
0595         lastPSH = &*ch;
0596         maxTOF = lastPSH->timeOfFlight();
0597       }
0598     }
0599   } else
0600     return result;  //return empty vector: no more hits on the sim track
0601   if (lastPSH == nullptr)
0602     return result;  //return empty vector: no more good hits on the sim track
0603   edm::LogVerbatim("CkfDebugger") << "CkfDebugger: corresponding SimHit is at gpos " << position(&*lastPSH);
0604 
0605   //take the simHits on the simTrack that are in the nextLayer (could be > 1 if overlap or matched)
0606   std::vector<PSimHit*> trackHits = idHitsMap[trajId];
0607   if (fabs((double)(trackHits.back()->detUnitId() - lastPSH->detUnitId())) < 1)
0608     return result;  //end of sim track
0609   std::vector<PSimHit*>::iterator currentIt = trackHits.end();
0610   for (std::vector<PSimHit*>::iterator it = trackHits.begin(); it != trackHits.end(); it++) {
0611     if (goodSimHit(**it) &&                                   //good hit
0612         (lastPSH->timeOfFlight() < (*it)->timeOfFlight()) &&  //greater tof
0613         //( fabs((double)((*it)->detUnitId()-(lastPSH->detUnitId()) ))>1) && //not components of the same matched hit
0614         ((det(lastPSH)->geographicalId().subdetId() != det(*it)->geographicalId().subdetId()) ||
0615          (layer(det(lastPSH)) != layer(det(*it))))  //change layer or detector(tib,tob,...)
0616     ) {
0617       edm::LogVerbatim("CkfDebugger") << "Next good PSimHit is at gpos " << position(*it);
0618       result.push_back(*it);
0619       currentIt = it;
0620       break;
0621     }
0622   }
0623   bool samelayer = true;
0624   if (currentIt != (trackHits.end() - 1) && currentIt != trackHits.end()) {
0625     for (std::vector<PSimHit*>::iterator nextIt = currentIt; (samelayer && nextIt != trackHits.end()); nextIt++) {
0626       if (goodSimHit(**nextIt)) {
0627         if ((det(*nextIt)->geographicalId().subdetId() == det(*currentIt)->geographicalId().subdetId()) &&
0628             (layer(det(*nextIt)) == layer(det(*currentIt)))) {
0629           result.push_back(*nextIt);
0630         } else
0631           samelayer = false;
0632       }
0633     }
0634   }
0635 
0636   return result;
0637 }
0638 
0639 bool CkfDebugger::goodSimHit(const PSimHit& sh) const {
0640   if (sh.pabs() > 0.9)
0641     return true;  // GeV, reject delta rays from association
0642   else
0643     return false;
0644 }
0645 
0646 bool CkfDebugger::associated(CTTRHp rechit, const PSimHit& pSimHit) const {
0647   LogTrace("CkfDebugger") << "now in associated";
0648 
0649   if (!rechit->isValid())
0650     return false;
0651   //   LogTrace("CkfDebugger") << "rec hit valid" ;
0652   const std::vector<PSimHit>& pSimHitVec = hitAssociator->associateHit(*rechit->hit());
0653   //   LogTrace("CkfDebugger") << "size=" << pSimHitVec.size() ;
0654   for (std::vector<PSimHit>::const_iterator shit = pSimHitVec.begin(); shit != pSimHitVec.end(); shit++) {
0655     //const GeomDetUnit* detUnit = theTrackerGeom->idToDetUnit( DetId(shit->detUnitId()));
0656     //         LogTrace("CkfDebugger") << "pSimHit.timeOfFlight()=" << pSimHit.timeOfFlight()
0657     //       << " pSimHit.pabs()=" << pSimHit.pabs() << " GP=" << position(&pSimHit);
0658     //         LogTrace("CkfDebugger") << "(*shit).timeOfFlight()=" << (*shit).timeOfFlight()
0659     //       << " (*shit).pabs()=" << (*shit).pabs() << " GP=" << detUnit->toGlobal( shit->localPosition());
0660     if ((fabs((*shit).timeOfFlight() - pSimHit.timeOfFlight()) < 1e-9) &&
0661         (fabs((*shit).pabs() - pSimHit.pabs()) < 1e-9))
0662       return true;
0663   }
0664   return false;
0665 }
0666 
0667 bool CkfDebugger::correctMeas(const TM& tm, const PSimHit* correctHit) const {
0668   LogTrace("CkfDebugger") << "now in correctMeas";
0669   const CTTRHp& recHit = tm.recHit();
0670   if (recHit->isValid())
0671     LogTrace("CkfDebugger") << "hit at position:" << recHit->globalPosition();
0672   TransientTrackingRecHit::RecHitContainer comp = recHit->transientHits();
0673   if (comp.empty()) {
0674     //     LogTrace("CkfDebugger") << "comp.empty()==true" ;
0675     return associated(recHit, *correctHit);
0676   } else {
0677     for (TransientTrackingRecHit::RecHitContainer::const_iterator ch = comp.begin(); ch != comp.end(); ++ch) {
0678       if (associated(recHit, *correctHit)) {
0679         // check if the other components are associated to the same trackId
0680         for (TransientTrackingRecHit::RecHitContainer::const_iterator ch2 = comp.begin(); ch2 != comp.end(); ++ch2) {
0681           if (ch2 == ch)
0682             continue;
0683           //////////
0684           //      LogTrace("CkfDebugger") << "correctHit->trackId()=" << correctHit->trackId() ;
0685           bool test = true;
0686           std::vector<SimHitIdpr> ids = hitAssociator->associateHitId(*(*ch2)->hit());
0687           for (std::vector<SimHitIdpr>::iterator j = ids.begin(); j != ids.end(); j++) {
0688             //      LogTrace("CkfDebugger") << "j=" <<j->first;
0689             if (correctHit->trackId() == j->first) {
0690               test = false;
0691               //          LogTrace("CkfDebugger") << correctHit->trackId()<< " " <<j->first;
0692             }
0693           }
0694           if (assocTrackId(*ch2) != ((int)(correctHit->trackId()))) {
0695             LogTrace("CkfDebugger") << "returning false 1"; /*return false;*/
0696           }  //fixme
0697           if (test) {
0698             //      LogTrace("CkfDebugger") << "returning false 2" ;
0699             return false;  // not all components from same simtrack
0700           }
0701           //      if (assocTrackId( **ch2) != ((int)( correctHit->trackId())) ) {
0702           //        return false; // not all components from same simtrack
0703           //      }
0704         }
0705         return true;  // if all components from same simtrack
0706       }
0707     }
0708     return false;
0709   }
0710 }
0711 
0712 //this checks only if there is the rechit on the det where the sim hit is
0713 pair<CTTRHp, double> CkfDebugger::analyseRecHitExistance(const PSimHit& sh, const TSOS& startingState) {
0714   LogTrace("CkfDebugger") << "now in analyseRecHitExistance";
0715 
0716 #if 0
0717   std::pair<CTTRHp, double> result;
0718   
0719   const MeasurementDet* simHitDet = theMeasurementTracker->idToDet( DetId( sh.detUnitId()));
0720   TSOS simHitState = TSOSFromSimHitFactory()(sh, *det(&sh), *theMagField);
0721   MeasurementDet::RecHitContainer recHits = simHitDet->recHits( simHitState);//take all hits from det
0722 
0723   //check if the hit is not present or is a problem of association
0724   TSOS firstDetState = theForwardPropagator->propagate( startingState, det(&sh)->surface());
0725   if (!firstDetState.isValid()) {
0726     edm::LogVerbatim("CkfDebugger") << "CkfDebugger: propagation failed from state " << startingState << " to first det surface " 
0727                     << position(&sh) ;
0728     propagation++;
0729     return std::pair<CTTRHp, double>((CTTRHp)(0),-1);
0730   }
0731 
0732   bool found = false;
0733   for ( MeasurementDet::RecHitContainer::const_iterator rh = recHits.begin(); rh != recHits.end(); rh++) {
0734     if ( associated( *rh, sh)) {
0735       found = true;
0736       result = std::pair<CTTRHp, double>(*rh,theChi2->estimate( firstDetState, **rh).second);
0737       edm::LogVerbatim("CkfDebugger") << "CkfDebugger: A RecHit associated to the correct Simhit exists at lpos " 
0738                       << (**rh).localPosition()
0739                       << " gpos " << (**rh).globalPosition()
0740                       << " layer " <<   layer((**rh).det())
0741                       << " subdet " << (**rh).det()->geographicalId().subdetId() 
0742                       << " Chi2 " << theChi2->estimate( firstDetState, **rh).second;
0743     }
0744   }
0745   if (!found) {
0746     edm::LogVerbatim("CkfDebugger") << "CkfDebugger: there is no RecHit associated to the correct SimHit." ;
0747     edm::LogVerbatim("CkfDebugger") << " There are " <<  recHits.size() << " RecHits in the simHit DetUnit" ;
0748     edm::LogVerbatim("CkfDebugger") << "SH GP=" << position(&sh) << " subdet=" << det(&sh)->geographicalId().subdetId() 
0749                     << " layer=" << layer(det(&sh)) ;
0750     int y=0;
0751     for (MeasurementDet::RecHitContainer::const_iterator rh = recHits.begin(); rh != recHits.end(); rh++)
0752       edm::LogVerbatim("CkfDebugger") << "RH#" << y++ << " GP=" << (**rh).globalPosition() << " subdet=" << (**rh).det()->geographicalId().subdetId() 
0753                       << " layer=" << layer((**rh).det()) ;
0754     for ( MeasurementDet::RecHitContainer::const_iterator rh = recHits.begin(); rh != recHits.end(); rh++) {
0755       edm::LogVerbatim("CkfDebugger") << "Non-associated RecHit at pos " << (**rh).localPosition() ;
0756     }
0757   }
0758 
0759   bool found2 = false;
0760   const PSimHit* sh2;
0761   StripSubdetector subdet( det(&sh)->geographicalId());
0762   if (!subdet.glued()) {
0763     edm::LogVerbatim("CkfDebugger") << "The DetUnit is not part of a GluedDet" ;
0764     if (found) {
0765       if (result.second>30){
0766     LogTrace("CkfDebugger") << "rh->parameters()=" << result.first->parameters() ;
0767     LogTrace("CkfDebugger") << "rh->parametersError()=" << result.first->parametersError() ;
0768     MeasurementExtractor me(firstDetState);
0769     AlgebraicVector r(result.first->parameters() - me.measuredParameters(*result.first));
0770     LogTrace("CkfDebugger") << "me.measuredParameters(**rh)=" << me.measuredParameters(*result.first) ;
0771     LogTrace("CkfDebugger") << "me.measuredError(**rh)=" << me.measuredError(*result.first) ;
0772     AlgebraicSymMatrix R(result.first->parametersError() + me.measuredError(*result.first));
0773     LogTrace("CkfDebugger") << "r=" << r ;
0774     LogTrace("CkfDebugger") << "R=" << R ;
0775     int ierr; 
0776     R.invert(ierr);
0777     LogTrace("CkfDebugger") << "R(-1)=" << R ;
0778     LogTrace("CkfDebugger") << "chi2=" << R.similarity(r) ;
0779       }
0780       return result;
0781     }
0782     else {
0783       simple_hit_not_found++;
0784       return std::pair<CTTRHp, double>((CTTRHp)(0),-8);//not found (is not a glued det)
0785     }
0786   } else {
0787     edm::LogVerbatim("CkfDebugger") << "The DetUnit is part of a GluedDet" ;
0788     DetId partnerDetId = DetId( subdet.partnerDetId());
0789 
0790     sh2 = pSimHit( sh.trackId(), partnerDetId);
0791     if (sh2 == 0) {
0792       edm::LogVerbatim("CkfDebugger") << "Partner DetUnit does not have a SimHit from the same track" ;
0793       if (found) {
0794     //projected rec hit
0795     TrackingRecHitProjector<ProjectedRecHit2D> proj;
0796     DetId gid = gluedId( subdet);
0797     const MeasurementDet* gluedDet = theMeasurementTracker->idToDet( gid);
0798     TSOS gluedTSOS = theForwardPropagator->propagate(startingState, gluedDet->geomDet().surface());
0799     CTTRHp projHit = proj.project( *result.first,gluedDet->geomDet(),gluedTSOS).get();
0800     //LogTrace("CkfDebugger") << proj.project( *result.first,gluedDet->geomDet(),gluedTSOS)->parameters() ;
0801     //LogTrace("CkfDebugger") << projHit->parametersError() ;
0802     double chi2 = theChi2->estimate(gluedTSOS, *proj.project( *result.first,gluedDet->geomDet(),gluedTSOS)).second;
0803     return std::pair<CTTRHp, double>(projHit,chi2);
0804       }
0805     }
0806     else {
0807       edm::LogVerbatim("CkfDebugger") << "Partner DetUnit has a good SimHit at gpos " << position(sh2) 
0808                       << " lpos " << sh2->localPosition() ;
0809       //}
0810     
0811       const MeasurementDet* partnerDet = theMeasurementTracker->idToDet( partnerDetId);
0812       if (partnerDet == 0) {
0813     edm::LogVerbatim("CkfDebugger") << "Partner measurementDet not found!!!" ;
0814     partner_det_not_fuond++;
0815     return std::pair<CTTRHp, double>((CTTRHp)(0),-3);
0816       }
0817       TSOS simHitState2 = TSOSFromSimHitFactory()(*sh2, *det(sh2), *theMagField);
0818       MeasurementDet::RecHitContainer recHits2 = partnerDet->recHits( simHitState2);
0819 
0820       TSOS secondDetState = theForwardPropagator->propagate( startingState, det(sh2)->surface());
0821       if (!secondDetState.isValid()) {
0822     edm::LogVerbatim("CkfDebugger") << "CkfDebugger: propagation failed from state " << startingState << " to second det surface " 
0823                     << position(sh2) ;
0824     propagation++;
0825     return std::pair<CTTRHp, double>((CTTRHp)(0),-1);
0826       }
0827 
0828       for ( MeasurementDet::RecHitContainer::const_iterator rh = recHits2.begin(); rh != recHits2.end(); rh++) {
0829     if ( associated( *rh, *sh2)) {
0830       found2 = true;
0831       edm::LogVerbatim("CkfDebugger") << "CkfDebugger: A RecHit associated to the correct Simhit exists at lpos " 
0832                       << (**rh).localPosition()
0833                       << " gpos " << (**rh).globalPosition()
0834                       << " Chi2 " << theChi2->estimate( secondDetState, **rh).second
0835         ;
0836     }
0837       }
0838       if (!found2) {
0839     edm::LogVerbatim("CkfDebugger") << "CkfDebugger: there is no RecHit associated to the correct SimHit." ;
0840     LogTrace("CkfDebugger") << " There are " <<  recHits.size() << " RecHits in the simHit DetUnit" ;
0841     for ( MeasurementDet::RecHitContainer::const_iterator rh = recHits.begin(); rh != recHits.end(); rh++) {
0842       LogTrace("CkfDebugger") << "Non-associated RecHit at pos " << (**rh).localPosition() ;
0843     }
0844       }
0845     }
0846   }
0847 
0848   MeasurementDet::RecHitContainer gluedHits;
0849   if (found && found2) {
0850     // look in the glued det
0851     DetId gid = gluedId( subdet);
0852     const MeasurementDet* gluedDet = theMeasurementTracker->idToDet( gid);
0853     if ( gluedDet == 0) {
0854       edm::LogVerbatim("CkfDebugger") << "CkfDebugger ERROR: glued MeasurementDet not found!" ;
0855       glued_det_not_fuond++;
0856       return std::pair<CTTRHp, double>((CTTRHp)(0),-4);
0857     }
0858 
0859     TSOS gluedDetState = theForwardPropagator->propagate( startingState, gluedDet->surface());
0860     if (!gluedDetState.isValid()) {
0861       edm::LogVerbatim("CkfDebugger") << "CkfDebugger: propagation failed from state " << startingState << " to det surface " 
0862                       << gluedDet->position() ;
0863       propagation++;
0864       return std::pair<CTTRHp, double>((CTTRHp)(0),-1);
0865     }
0866 
0867     gluedHits = gluedDet->recHits( gluedDetState);
0868     edm::LogVerbatim("CkfDebugger") << "CkfDebugger: the GluedDet returned " << gluedHits.size() << " hits" ;
0869     if (gluedHits.size()==0){
0870       edm::LogVerbatim("CkfDebugger") << "Found and associated mono and stereo recHits but not matched!!!" ;
0871       matched_not_found++;
0872       return std::pair<CTTRHp, double>((CTTRHp)(0),-5);
0873     } 
0874     bool found3 = false;
0875     for ( MeasurementDet::RecHitContainer::const_iterator rh = gluedHits.begin(); rh != gluedHits.end(); rh++) {
0876       if ( associated( *rh, sh) && associated( *rh, *sh2)) {
0877     double chi2 = theChi2->estimate(gluedDetState, **rh).second;
0878     edm::LogVerbatim("CkfDebugger") << "Matched hit at lpos " << (**rh).localPosition()
0879                     << " gpos " << (**rh).globalPosition()
0880                     << " has Chi2 " << chi2
0881       ;
0882     result = std::pair<CTTRHp, double>(&**rh,chi2);
0883     found3 = true;
0884     if (chi2>30){
0885       LogTrace("CkfDebugger") << "rh->parameters()=" << (*rh)->parameters() ;
0886       LogTrace("CkfDebugger") << "rh->parametersError()=" << (*rh)->parametersError() ;
0887       MeasurementExtractor me(gluedDetState);
0888       AlgebraicVector r((*rh)->parameters() - me.measuredParameters(**rh));
0889       LogTrace("CkfDebugger") << "me.measuredParameters(**rh)=" << me.measuredParameters(**rh) ;
0890       LogTrace("CkfDebugger") << "me.measuredError(**rh)=" << me.measuredError(**rh) ;
0891       AlgebraicSymMatrix R((*rh)->parametersError() + me.measuredError(**rh));
0892       LogTrace("CkfDebugger") << "r=" << r ;
0893       LogTrace("CkfDebugger") << "R=" << R ;
0894       int ierr; 
0895       R.invert(ierr);
0896       LogTrace("CkfDebugger") << "R(-1)=" << R ;
0897       LogTrace("CkfDebugger") << "chi2=" << R.similarity(r) ;
0898     }
0899     break;
0900       }
0901     }
0902     if (found3) return result;
0903     else {
0904       edm::LogVerbatim("CkfDebugger") << "Found and associated mono and stereo recHits. Matched found but not associated!!!" ;
0905       matched_not_associated++;
0906       return std::pair<CTTRHp, double>((CTTRHp)(0),-6);
0907     }
0908   }
0909   else if ( (found && !found2) || (!found && found2) ) {
0910     edm::LogVerbatim("CkfDebugger") << "Only one component is found" ;
0911     only_one_component++;
0912     return std::pair<CTTRHp, double>((CTTRHp)(0),-7);
0913   }
0914   else {
0915     edm::LogVerbatim("CkfDebugger") << "No component is found" ;
0916     no_component++;
0917     return std::pair<CTTRHp, double>((CTTRHp)(0),-2);
0918   }
0919   other++;
0920 #endif
0921   return std::pair<CTTRHp, double>((CTTRHp)(nullptr), 0);  //other
0922 }
0923 
0924 const PSimHit* CkfDebugger::pSimHit(unsigned int tkId, DetId detId) {
0925   for (std::vector<PSimHit*>::iterator shi = idHitsMap[tkId].begin(); shi != idHitsMap[tkId].end(); ++shi) {
0926     if ((*shi)->detUnitId() == detId.rawId() &&
0927         //(shi)->trackId() == tkId &&
0928         goodSimHit(**shi)) {
0929       return (*shi);
0930     }
0931   }
0932   return nullptr;
0933 }
0934 
0935 int CkfDebugger::analyseRecHitNotFound(const Trajectory& traj, CTTRHp correctRecHit) {
0936   unsigned int correctDetId = correctRecHit->det()->geographicalId().rawId();
0937   int correctLayId = layer(correctRecHit->det());
0938   LogTrace("CkfDebugger") << "correct layer id=" << correctLayId;
0939 
0940   TSOS currentState(traj.lastMeasurement().updatedState());
0941   std::vector<const DetLayer*> nl =
0942       theNavSchool->nextLayers(*traj.lastLayer(), *currentState.freeState(), traj.direction());
0943   if (nl.empty()) {
0944     edm::LogVerbatim("CkfDebugger") << "no compatible layers";
0945     no_layer++;
0946     return 2;
0947   }
0948 
0949   TkLayerLess lless;  //FIXME - was lless(traj.direction())
0950   const DetLayer* detLayer = nullptr;
0951   bool navLayerAfter = false;
0952   bool test = false;
0953   for (std::vector<const DetLayer*>::iterator il = nl.begin(); il != nl.end(); il++) {
0954     if (dynamic_cast<const BarrelDetLayer*>(*il)) {
0955       const BarrelDetLayer* pbl = dynamic_cast<const BarrelDetLayer*>(*il);
0956       LogTrace("CkfDebugger") << "pbl->specificSurface().bounds().length()="
0957                               << pbl->specificSurface().bounds().length();
0958       LogTrace("CkfDebugger") << "pbl->specificSurface().bounds().width()=" << pbl->specificSurface().bounds().width();
0959     }
0960     int layId = layer(((*(*il)->basicComponents().begin())));
0961     LogTrace("CkfDebugger") << " subdet=" << (*(*il)->basicComponents().begin())->geographicalId().subdetId()
0962                             << "layer id=" << layId;
0963     if (layId == correctLayId) {
0964       test = true;
0965       detLayer = &**il;
0966       break;
0967     }
0968     if (lless(*il, theGeomSearchTracker->detLayer(correctRecHit->det()->geographicalId())))
0969       navLayerAfter = true;  //it is enough that only one layer is after the correct one?
0970   }
0971 
0972   if (test) {
0973     edm::LogVerbatim("CkfDebugger") << "correct layer taken into account. layer id: " << correctLayId;
0974   } else if (navLayerAfter) {
0975     edm::LogVerbatim("CkfDebugger") << "SimHit layer after the layers returned by Navigation.";
0976     edm::LogVerbatim("CkfDebugger") << "Probably a missing SimHit.";
0977     edm::LogVerbatim("CkfDebugger") << "check: " << (correctRecHit->det()->geographicalId().subdetId()) << " "
0978                                     << (layer(correctRecHit->det()));
0979     dump6[pair<int, int>((correctRecHit->det()->geographicalId().subdetId() - 1), (layer(correctRecHit->det())) - 1)]++;
0980     no_sim_hit++;
0981     return 16;
0982   } else {
0983     edm::LogVerbatim("CkfDebugger") << "correct layer NOT taken into account. correct layer id: " << correctLayId;
0984     layer_not_found++;
0985     return 3;
0986   }
0987 
0988   typedef DetLayer::DetWithState DetWithState;
0989   std::vector<DetWithState> compatDets = detLayer->compatibleDets(currentState, *theForwardPropagator, *theChi2);
0990   //   LogTrace("CkfDebugger") << "DEBUGGER" ;
0991   //   LogTrace("CkfDebugger") << "runned compatDets." ;
0992   //   LogTrace("CkfDebugger") << "started from the following TSOS:" ;
0993   //   LogTrace("CkfDebugger") << currentState ;
0994   //   LogTrace("CkfDebugger") << "number of dets found=" << compatDets.size() ;
0995   //   for (std::vector<DetWithState>::iterator det=compatDets.begin();det!=compatDets.end();det++){
0996   //     unsigned int detId = det->first->geographicalId().rawId();
0997   //     LogTrace("CkfDebugger") << "detId=" << detId ;
0998   //   }
0999   bool test2 = false;
1000   for (std::vector<DetWithState>::iterator det = compatDets.begin(); det != compatDets.end(); det++) {
1001     unsigned int detId = det->first->geographicalId().rawId();
1002     //     LogTrace("CkfDebugger") << "detId=" << detId
1003     //   << "\ncorrectRecHit->det()->geographicalId().rawId()=" << correctRecHit->det()->geographicalId().rawId()
1004     //   << "\ngluedId(correctRecHit->det()->geographicalId()).rawId()=" << gluedId(correctRecHit->det()->geographicalId()).rawId()
1005     //   ;
1006     if (detId == gluedId(correctRecHit->det()->geographicalId()).rawId()) {
1007       test2 = true;
1008       break;
1009     }
1010   }
1011 
1012   if (test2) {
1013     edm::LogVerbatim("CkfDebugger") << "correct det taken into account. correctDetId is: " << correctDetId
1014                                     << ". please check chi2.";
1015     return 5;
1016   } else {
1017     edm::LogVerbatim("CkfDebugger") << "correct det NOT taken into account. correctDetId: " << correctDetId;
1018     det_not_found++;
1019     return 4;
1020   }
1021 }
1022 
1023 double CkfDebugger::testSeed(CTTRHp recHit1, CTTRHp recHit2, TSOS state) {
1024   //edm::LogVerbatim("CkfDebugger") << "CkfDebugger::testSeed";
1025   //test Deltas
1026   const std::vector<PSimHit>& pSimHitVec1 = hitAssociator->associateHit(*recHit1->hit());
1027   const std::vector<PSimHit>& pSimHitVec2 = hitAssociator->associateHit(*recHit2->hit());
1028 
1029   if (pSimHitVec1.empty() || pSimHitVec2.empty() || hasDelta(&(*pSimHitVec1.begin())) ||
1030       hasDelta(&(*pSimHitVec2.begin()))) {
1031     edm::LogVerbatim("CkfDebugger") << "Seed has delta or problems";
1032     return -1;
1033   }
1034 
1035   //   LogTrace("CkfDebugger") << "state=\n" << state ;
1036   //   double stlp1 = state.localParameters().vector()[0];
1037   //   double stlp2 = state.localParameters().vector()[1];
1038   //   double stlp3 = state.localParameters().vector()[2];
1039   //   double stlp4 = state.localParameters().vector()[3];
1040   //   double stlp5 = state.localParameters().vector()[4];
1041 
1042   if (!pSimHitVec2.empty()) {
1043     const PSimHit& simHit = *pSimHitVec2.begin();
1044 
1045     double shlp1 = -1 / simHit.momentumAtEntry().mag();
1046     double shlp2 = simHit.momentumAtEntry().x() / simHit.momentumAtEntry().z();
1047     double shlp3 = simHit.momentumAtEntry().y() / simHit.momentumAtEntry().z();
1048     double shlp4 = simHit.localPosition().x();
1049     double shlp5 = simHit.localPosition().y();
1050     AlgebraicVector5 v;
1051     v[0] = shlp1;
1052     v[1] = shlp2;
1053     v[2] = shlp3;
1054     v[3] = shlp4;
1055     v[4] = shlp5;
1056 
1057     //     LogTrace("CkfDebugger") << "simHit.localPosition()=" << simHit.localPosition() ;
1058     //     LogTrace("CkfDebugger") << "simHit.momentumAtEntry()=" << simHit.momentumAtEntry() ;
1059     //     LogTrace("CkfDebugger") << "recHit2->localPosition()=" << recHit2->localPosition() ;
1060     //     LogTrace("CkfDebugger") << "recHit2->localPositionError()=" << recHit2->localPositionError() ;
1061     //     LogTrace("CkfDebugger") << "state.localPosition()=" << state.localPosition() ;
1062     //     LogTrace("CkfDebugger") << "state.localError().positionError()=" << state.localError().positionError() ;
1063 
1064     //     LogTrace("CkfDebugger") << "pullx(sh-rh)=" << (simHit.localPosition().x()-recHit2->localPosition().x())/sqrt(recHit2->localPositionError().xx()) ;
1065     //     LogTrace("CkfDebugger") << "pullx(sh-st)=" << (simHit.localPosition().x()-state.localPosition().x())/sqrt(state.localError().positionError().xx()) ;
1066     //     LogTrace("CkfDebugger") << "pullx(st-rh)=" << (state.localPosition().x()-recHit2->localPosition().x())/
1067     //       sqrt(recHit2->localPositionError().xx()+state.localError().positionError().xx()) ;
1068 
1069     //     LogTrace("CkfDebugger") << "local parameters" ;
1070     //     LogTrace("CkfDebugger") << left;
1071     //     LogTrace("CkfDebugger") << setw(15) << stlp1 << setw(15) << shlp1 << setw(15) << sqrt(state.localError().matrix()[0][0])
1072     //   << setw(15) << (stlp1-shlp1)/stlp1 << setw(15) << (stlp1-shlp1)/sqrt(state.localError().matrix()[0][0]) ;
1073     //     LogTrace("CkfDebugger") << setw(15) << stlp2 << setw(15) << shlp2 << setw(15) << sqrt(state.localError().matrix()[1][1])
1074     //   << setw(15) << (stlp2-shlp2)/stlp2 << setw(15) << (stlp2-shlp2)/sqrt(state.localError().matrix()[1][1]) ;
1075     //     LogTrace("CkfDebugger") << setw(15) << stlp3 << setw(15) << shlp3 << setw(15) << sqrt(state.localError().matrix()[2][2])
1076     //       << setw(15) << (stlp3-shlp3)/stlp3 << setw(15) << (stlp3-shlp3)/sqrt(state.localError().matrix()[2][2]) ;
1077     //     LogTrace("CkfDebugger") << setw(15) << stlp4 << setw(15) << shlp4 << setw(15) << sqrt(state.localError().matrix()[3][3])
1078     //       << setw(15) << (stlp4-shlp4)/stlp4 << setw(15) << (stlp4-shlp4)/sqrt(state.localError().matrix()[3][3]) ;
1079     //     LogTrace("CkfDebugger") << setw(15) << stlp5 << setw(15) << shlp5 << setw(15) << sqrt(state.localError().matrix()[4][4]) <<
1080     //       setw(15) << (stlp5-shlp5)/stlp5 << setw(15) << (stlp5-shlp5)/sqrt(state.localError().matrix()[4][4]) ;
1081 
1082     AlgebraicSymMatrix55 R = state.localError().matrix();
1083     R.Invert();
1084     double chi2 = ROOT::Math::Similarity(v - state.localParameters().vector(), R);
1085     LogTrace("CkfDebugger") << "chi2=" << chi2;
1086     return chi2;
1087   }
1088 
1089   return 0;  //fixme
1090 }
1091 
1092 CkfDebugger::~CkfDebugger() {
1093   for (int it = 0; it != ((int)(dump.size())); it++)
1094     edm::LogVerbatim("CkfDebugger") << "dump " << it << " " << dump[it];
1095 
1096   edm::LogVerbatim("CkfDebugger");
1097   edm::LogVerbatim("CkfDebugger") << "seedWithDelta=" << ((double)seedWithDelta / totSeeds);
1098   edm::LogVerbatim("CkfDebugger") << "problems=" << ((double)problems / totSeeds);
1099   edm::LogVerbatim("CkfDebugger") << "no_sim_hit=" << ((double)no_sim_hit / totSeeds);
1100   edm::LogVerbatim("CkfDebugger") << "no_layer=" << ((double)no_layer / totSeeds);
1101   edm::LogVerbatim("CkfDebugger") << "layer_not_found=" << ((double)layer_not_found / totSeeds);
1102   edm::LogVerbatim("CkfDebugger") << "det_not_found=" << ((double)det_not_found / totSeeds);
1103   edm::LogVerbatim("CkfDebugger") << "chi2gt30=" << ((double)chi2gt30 / totSeeds);
1104   edm::LogVerbatim("CkfDebugger") << "chi2gt30deltaSeed=" << ((double)chi2gt30deltaSeed / totSeeds);
1105   edm::LogVerbatim("CkfDebugger") << "chi2gt30delta=" << ((double)chi2gt30delta / totSeeds);
1106   edm::LogVerbatim("CkfDebugger") << "chi2ls30=" << ((double)chi2ls30 / totSeeds);
1107   edm::LogVerbatim("CkfDebugger") << "simple_hit_not_found=" << ((double)simple_hit_not_found / totSeeds);
1108   edm::LogVerbatim("CkfDebugger") << "no_component=" << ((double)no_component / totSeeds);
1109   edm::LogVerbatim("CkfDebugger") << "only_one_component=" << ((double)only_one_component / totSeeds);
1110   edm::LogVerbatim("CkfDebugger") << "matched_not_found=" << ((double)matched_not_found / totSeeds);
1111   edm::LogVerbatim("CkfDebugger") << "matched_not_associated=" << ((double)matched_not_associated / totSeeds);
1112   edm::LogVerbatim("CkfDebugger") << "partner_det_not_fuond=" << ((double)partner_det_not_fuond / totSeeds);
1113   edm::LogVerbatim("CkfDebugger") << "glued_det_not_fuond=" << ((double)glued_det_not_fuond / totSeeds);
1114   edm::LogVerbatim("CkfDebugger") << "propagation=" << ((double)propagation / totSeeds);
1115   edm::LogVerbatim("CkfDebugger") << "other=" << ((double)other / totSeeds);
1116   edm::LogVerbatim("CkfDebugger") << "totchi2gt30=" << ((double)totchi2gt30 / totSeeds);
1117   edm::LogVerbatim("CkfDebugger") << "totSeeds=" << totSeeds;
1118   edm::LogVerbatim("CkfDebugger");
1119 
1120   edm::LogVerbatim("CkfDebugger") << "layer navigation problems:";
1121   for (int i = 0; i != 6; i++)
1122     for (int j = 0; j != 9; j++) {
1123       if (i == 0 && j > 2)
1124         break;
1125       if (i == 1 && j > 1)
1126         break;
1127       if (i == 2 && j > 3)
1128         break;
1129       if (i == 3 && j > 2)
1130         break;
1131       if (i == 4 && j > 5)
1132         break;
1133       if (i == 5 && j > 8)
1134         break;
1135       edm::LogVerbatim("CkfDebugger") << "det=" << i + 1 << " lay=" << j + 1 << " " << dump2[pair<int, int>(i, j)];
1136     }
1137   edm::LogVerbatim("CkfDebugger") << "\nlayer with hit having chi2>30:";
1138   for (int i = 0; i != 6; i++)
1139     for (int j = 0; j != 9; j++) {
1140       if (i == 0 && j > 2)
1141         break;
1142       if (i == 1 && j > 1)
1143         break;
1144       if (i == 2 && j > 3)
1145         break;
1146       if (i == 3 && j > 2)
1147         break;
1148       if (i == 4 && j > 5)
1149         break;
1150       if (i == 5 && j > 8)
1151         break;
1152       edm::LogVerbatim("CkfDebugger") << "det=" << i + 1 << " lay=" << j + 1 << " " << dump3[pair<int, int>(i, j)];
1153     }
1154   edm::LogVerbatim("CkfDebugger") << "\nlayer with hit having chi2>30 for delta rays:";
1155   for (int i = 0; i != 6; i++)
1156     for (int j = 0; j != 9; j++) {
1157       if (i == 0 && j > 2)
1158         break;
1159       if (i == 1 && j > 1)
1160         break;
1161       if (i == 2 && j > 3)
1162         break;
1163       if (i == 3 && j > 2)
1164         break;
1165       if (i == 4 && j > 5)
1166         break;
1167       if (i == 5 && j > 8)
1168         break;
1169       edm::LogVerbatim("CkfDebugger") << "det=" << i + 1 << " lay=" << j + 1 << " " << dump5[pair<int, int>(i, j)];
1170     }
1171   edm::LogVerbatim("CkfDebugger") << "\nlayer with det not found:";
1172   for (int i = 0; i != 6; i++)
1173     for (int j = 0; j != 9; j++) {
1174       if (i == 0 && j > 2)
1175         break;
1176       if (i == 1 && j > 1)
1177         break;
1178       if (i == 2 && j > 3)
1179         break;
1180       if (i == 3 && j > 2)
1181         break;
1182       if (i == 4 && j > 5)
1183         break;
1184       if (i == 5 && j > 8)
1185         break;
1186       edm::LogVerbatim("CkfDebugger") << "det=" << i + 1 << " lay=" << j + 1 << " " << dump4[pair<int, int>(i, j)];
1187     }
1188   edm::LogVerbatim("CkfDebugger") << "\nlayer with correct RecHit after missing Sim Hit:";
1189   for (int i = 0; i != 6; i++)
1190     for (int j = 0; j != 9; j++) {
1191       if (i == 0 && j > 2)
1192         break;
1193       if (i == 1 && j > 1)
1194         break;
1195       if (i == 2 && j > 3)
1196         break;
1197       if (i == 3 && j > 2)
1198         break;
1199       if (i == 4 && j > 5)
1200         break;
1201       if (i == 5 && j > 8)
1202         break;
1203       edm::LogVerbatim("CkfDebugger") << "det=" << i + 1 << " lay=" << j + 1 << " " << dump6[pair<int, int>(i, j)];
1204     }
1205   hchi2seedAll->Write();
1206   hchi2seedProb->Write();
1207   std::stringstream title;
1208   for (int i = 0; i != 6; i++)
1209     for (int j = 0; j != 9; j++) {
1210       if (i == 0 && j > 2)
1211         break;
1212       if (i == 1 && j > 1)
1213         break;
1214       if (i == 2 && j > 3)
1215         break;
1216       if (i == 3 && j > 2)
1217         break;
1218       if (i == 4 && j > 5)
1219         break;
1220       if (i == 5 && j > 8)
1221         break;
1222       title.str("");
1223       title << "pullX_" << i + 1 << "-" << j + 1 << "_sh-rh";
1224       hPullX_shrh[title.str()]->Write();
1225       title.str("");
1226       title << "pullY_" << i + 1 << "-" << j + 1 << "_sh-rh";
1227       hPullY_shrh[title.str()]->Write();
1228       title.str("");
1229       title << "pullX_" << i + 1 << "-" << j + 1 << "_sh-st";
1230       hPullX_shst[title.str()]->Write();
1231       title.str("");
1232       title << "pullY_" << i + 1 << "-" << j + 1 << "_sh-st";
1233       hPullY_shst[title.str()]->Write();
1234       title.str("");
1235       title << "pullX_" << i + 1 << "-" << j + 1 << "_st-rh";
1236       hPullX_strh[title.str()]->Write();
1237       title.str("");
1238       title << "pullY_" << i + 1 << "-" << j + 1 << "_st-rh";
1239       hPullY_strh[title.str()]->Write();
1240       title.str("");
1241       title << "PullGP_X_" << i + 1 << "-" << j + 1 << "_sh-st";
1242       hPullGP_X_shst[title.str()]->Write();
1243       title.str("");
1244       title << "PullGP_Y_" << i + 1 << "-" << j + 1 << "_sh-st";
1245       hPullGP_Y_shst[title.str()]->Write();
1246       title.str("");
1247       title << "PullGP_Z_" << i + 1 << "-" << j + 1 << "_sh-st";
1248       hPullGP_Z_shst[title.str()]->Write();
1249       if (((i == 2 || i == 4) && (j == 0 || j == 1)) || (i == 3 || i == 5)) {
1250         title.str("");
1251         title << "pullM_" << i + 1 << "-" << j + 1 << "_sh-rh";
1252         hPullM_shrh[title.str()]->Write();
1253         title.str("");
1254         title << "pullS_" << i + 1 << "-" << j + 1 << "_sh-rh";
1255         hPullS_shrh[title.str()]->Write();
1256         title.str("");
1257         title << "pullM_" << i + 1 << "-" << j + 1 << "_sh-st";
1258         hPullM_shst[title.str()]->Write();
1259         title.str("");
1260         title << "pullS_" << i + 1 << "-" << j + 1 << "_sh-st";
1261         hPullS_shst[title.str()]->Write();
1262         title.str("");
1263         title << "pullM_" << i + 1 << "-" << j + 1 << "_st-rh";
1264         hPullM_strh[title.str()]->Write();
1265         title.str("");
1266         title << "pullS_" << i + 1 << "-" << j + 1 << "_st-rh";
1267         hPullS_strh[title.str()]->Write();
1268       }
1269     }
1270   hPullGPXvsGPX_shst->Write();
1271   hPullGPXvsGPY_shst->Write();
1272   hPullGPXvsGPZ_shst->Write();
1273   hPullGPXvsGPr_shst->Write();
1274   hPullGPXvsGPeta_shst->Write();
1275   hPullGPXvsGPphi_shst->Write();
1276 
1277   //file->Write();
1278   file->Close();
1279 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!