00001 #include "RecoTracker/DebugTools/interface/CkfDebugger.h"
00002 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
00003 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00004 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
00005 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00006 #include "Geometry/CommonTopologies/interface/Topology.h"
00007 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementVector.h"
00008 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
00009 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00010 #include "RecoTracker/DebugTools/interface/TSOSFromSimHitFactory.h"
00011 #include "TrackingTools/MeasurementDet/interface/MeasurementDet.h"
00012 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
00013 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
00014 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00015
00016 #include "RecoTracker/DebugTools/interface/GluedDetFromDetUnit.h"
00017 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00018
00019 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00020 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
00021 #include "RecoTracker/TkNavigation/interface/TkLayerLess.h"
00022 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
00023 #include "RecoTracker/TkDetLayers/interface/PixelBarrelLayer.h"
00024
00025 #include "TrackingTools/TransientTrackingRecHit/interface/TrackingRecHitProjector.h"
00026 #include "RecoTracker/TransientTrackingRecHit/interface/ProjectedRecHit2D.h"
00027
00028 #include <iostream>
00029 #include <iomanip>
00030 #include <sstream>
00031
00032 using namespace std;
00033
00034 CkfDebugger::CkfDebugger( edm::EventSetup const & es ):totSeeds(0)
00035 {
00036 file = new TFile("out.root","recreate");
00037 hchi2seedAll = new TH1F("hchi2seedAll","hchi2seedAll",2000,0,200);
00038 hchi2seedProb = new TH1F("hchi2seedProb","hchi2seedProb",2000,0,200);
00039
00040 edm::ESHandle<TrackerGeometry> tracker;
00041 es.get<TrackerDigiGeometryRecord>().get(tracker);
00042 theTrackerGeom = &(*tracker);
00043
00044 edm::ESHandle<MagneticField> theField;
00045 es.get<IdealMagneticFieldRecord>().get(theField);
00046 theMagField = &(*theField);
00047
00048 for (int i=0; i!=17; i++){
00049 dump.push_back(0);
00050 }
00051
00052 std::stringstream title;
00053 for (int i=0; i!=6; i++)
00054 for (int j=0; j!=9; j++){
00055 if (i==0 && j>2) break;
00056 if (i==1 && j>1) break;
00057 if (i==2 && j>3) break;
00058 if (i==3 && j>2) break;
00059 if (i==4 && j>5) break;
00060 if (i==5 && j>8) break;
00061 dump2[pair<int,int>(i,j)]=0;
00062 dump3[pair<int,int>(i,j)]=0;
00063 dump4[pair<int,int>(i,j)]=0;
00064 dump5[pair<int,int>(i,j)]=0;
00065 dump6[pair<int,int>(i,j)]=0;
00066 title.str("");
00067 title << "pullX_" << i+1 << "-" << j+1 << "_sh-rh";
00068 hPullX_shrh[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00069 title.str("");
00070 title << "pullY_" << i+1 << "-" << j+1 << "_sh-rh";
00071 hPullY_shrh[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00072 title.str("");
00073 title << "pullX_" << i+1 << "-" << j+1 << "_sh-st";
00074 hPullX_shst[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00075 title.str("");
00076 title << "pullY_" << i+1 << "-" << j+1 << "_sh-st";
00077 hPullY_shst[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00078 title.str("");
00079 title << "pullX_" << i+1 << "-" << j+1 << "_st-rh";
00080 hPullX_strh[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00081 title.str("");
00082 title << "pullY_" << i+1 << "-" << j+1 << "_st-rh";
00083 hPullY_strh[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00084 title.str("");
00085 title << "PullGP_X_" << i+1 << "-" << j+1 << "_sh-st";
00086 hPullGP_X_shst[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00087 title.str("");
00088 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_sh-st";
00089 hPullGP_Y_shst[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00090 title.str("");
00091 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_sh-st";
00092 hPullGP_Z_shst[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00093 if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
00094 title.str("");
00095 title << "pullM_" << i+1 << "-" << j+1 << "_sh-rh";
00096 hPullM_shrh[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00097 title.str("");
00098 title << "pullS_" << i+1 << "-" << j+1 << "_sh-rh";
00099 hPullS_shrh[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00100 title.str("");
00101 title << "pullM_" << i+1 << "-" << j+1 << "_sh-st";
00102 hPullM_shst[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00103 title.str("");
00104 title << "pullS_" << i+1 << "-" << j+1 << "_sh-st";
00105 hPullS_shst[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00106 title.str("");
00107 title << "pullM_" << i+1 << "-" << j+1 << "_st-rh";
00108 hPullM_strh[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00109 title.str("");
00110 title << "pullS_" << i+1 << "-" << j+1 << "_st-rh";
00111 hPullS_strh[title.str()] = new TH1F(title.str().c_str(),title.str().c_str(),1000,-50,50);
00112 }
00113 }
00114
00115 hPullGPXvsGPX_shst = new TH2F("PullGPXvsGPX_shst","PullGPXvsGPX_shst",1000,-50,50,100,-50,50);
00116 hPullGPXvsGPY_shst = new TH2F("PullGPXvsGPY_shst","PullGPXvsGPY_shst",1000,-50,50,100,-50,50);
00117 hPullGPXvsGPZ_shst = new TH2F("PullGPXvsGPZ_shst","PullGPXvsGPZ_shst",1000,-50,50,200,-100,100);
00118 hPullGPXvsGPr_shst = new TH2F("PullGPXvsGPr_shst","PullGPXvsGPr_shst",1000,-50,50,300,-150,150);
00119 hPullGPXvsGPeta_shst = new TH2F("PullGPXvsGPeta_shst","PullGPXvsGPeta_shst",1000,-50,50,50,-2.5,2.5);
00120 hPullGPXvsGPphi_shst = new TH2F("PullGPXvsGPphi_shst","PullGPXvsGPphi_shst",1000,-50,50,63,0,6.3);
00121
00122 seedWithDelta=0;
00123 problems=0;
00124 no_sim_hit=0;
00125 no_layer=0;
00126 layer_not_found=0;
00127 det_not_found=0;
00128 chi2gt30=0;
00129 chi2gt30delta=0;
00130 chi2gt30deltaSeed=0;
00131 chi2ls30=0;
00132 simple_hit_not_found=0;
00133 no_component=0;
00134 only_one_component=0;
00135 matched_not_found=0;
00136 matched_not_associated=0;
00137 partner_det_not_fuond=0;
00138 glued_det_not_fuond=0;
00139 propagation=0;
00140 other=0;
00141 totchi2gt30=0;
00142 }
00143
00144 void CkfDebugger::printSimHits( const edm::Event& iEvent)
00145 {
00146 edm::LogVerbatim("CkfDebugger") << "\nEVENT #" << iEvent.id();
00147
00148 hitAssociator = new TrackerHitAssociator(iEvent);
00149
00150 std::map<unsigned int, std::vector<PSimHit> >& theHitsMap = hitAssociator->SimHitMap;
00151 idHitsMap.clear();
00152
00153 for (std::map<unsigned int, std::vector<PSimHit> >::iterator it=theHitsMap.begin();
00154 it!=theHitsMap.end();it++){
00155 for (std::vector<PSimHit>::iterator isim = it->second.begin();
00156 isim != it->second.end(); ++isim){
00157 idHitsMap[isim->trackId()].push_back(&*isim);
00158 }
00159 }
00160
00161 for (std::map<unsigned int,std::vector<PSimHit*> >::iterator it=idHitsMap.begin();
00162 it!=idHitsMap.end();it++){
00163 sort(it->second.begin(),it->second.end(),less_mag());
00164 for (std::vector<PSimHit*>::iterator isim = it->second.begin();
00165 isim != it->second.end(); ++isim){
00166 const GeomDetUnit* detUnit = theTrackerGeom->idToDetUnit( DetId((*isim)->detUnitId()));
00167 dumpSimHit(SimHit( (*isim), detUnit));
00168 }
00169 }
00170
00171 }
00172
00173 void CkfDebugger::dumpSimHit(const SimHit& hit) const
00174 {
00175 GlobalPoint pos = hit.globalPosition();
00176 edm::LogVerbatim("CkfDebugger") << "SimHit pos" << pos
00177 << " r=" << pos.perp() << " phi=" << pos.phi()
00178 << " trackId=" << hit.trackId()
00179 << " particleType=" << hit.particleType()
00180 << " pabs=" << hit.pabs()
00181 << " processType=" << hit. processType();
00182 }
00183
00184
00185 bool CkfDebugger::analyseCompatibleMeasurements(const Trajectory& traj,
00186 const std::vector<TrajectoryMeasurement>& meas,
00187 const MeasurementTracker* aMeasurementTracker,
00188 const Propagator* propagator,
00189 const Chi2MeasurementEstimatorBase* estimator,
00190 const TransientTrackingRecHitBuilder* aTTRHBuilder)
00191 {
00192 LogTrace("CkfDebugger") << "\nnow in analyseCompatibleMeasurements" ;
00193 LogTrace("CkfDebugger") << "number of input hits:" << meas.size() ;
00194 for(std::vector<TrajectoryMeasurement>::const_iterator tmpIt=meas.begin();tmpIt!=meas.end();tmpIt++){
00195 if (tmpIt->recHit()->isValid()) LogTrace("CkfDebugger") << "valid hit at position:" << tmpIt->recHit()->globalPosition() ;
00196 }
00197 theForwardPropagator = propagator;
00198 theChi2 = estimator;
00199 theMeasurementTracker = aMeasurementTracker;
00200 theGeomSearchTracker = theMeasurementTracker->geometricSearchTracker();
00201 theTTRHBuilder = aTTRHBuilder;
00202 unsigned int trajId = 0;
00203 if ( !correctTrajectory(traj, trajId)) {
00204 LogTrace("CkfDebugger") << "trajectory not correct" ;
00205 return true;
00206 }
00207 LogTrace("CkfDebugger") << "correct trajectory" ;
00208
00209 if (traj.measurements().size() == 2){
00210 if ( testSeed(traj.firstMeasurement().recHit(),traj.lastMeasurement().recHit(),traj.lastMeasurement().updatedState()) == -1 ) {
00211 LogTrace("CkfDebugger") << "Seed has delta" ;
00212 seedWithDelta++;
00213 return false;
00214 }
00215 }
00216
00217
00218
00219 std::vector<const PSimHit*> correctHits = nextCorrectHits(traj, trajId);
00220 if ( correctHits.size() == 0) return true;
00221
00222 for (std::vector<const PSimHit*>::iterator corHit=correctHits.begin();corHit!=correctHits.end();corHit++){
00223 for (std::vector<TM>::const_iterator i=meas.begin(); i!=meas.end(); i++) {
00224 if (correctMeas( *i, *corHit)) {
00225 LogTrace("CkfDebugger") << "Correct hit found at position " << i-meas.begin() ;
00226 return true;
00227 }
00228 }
00229 }
00230
00231
00232
00233 const PSimHit* correctHit = *(correctHits.begin());
00234
00235
00236 edm::LogVerbatim("CkfDebugger") << std::endl << "CkfDebugger: problem found: correct hit not found by findCompatibleMeasurements" ;
00237 edm::LogVerbatim("CkfDebugger") << "The correct hit position is " << position(correctHit) << " lp " << correctHit->localPosition() ;
00238 edm::LogVerbatim("CkfDebugger") << "The size of the meas vector is " << meas.size() ;
00239 dump[0]++;problems++;
00240
00241 for (std::vector<TM>::const_iterator i=meas.begin(); i!=meas.end(); i++) {
00242 edm::LogVerbatim("CkfDebugger") << "Is the hit valid? " << i->recHit()->isValid() ;
00243 if (i->recHit()->isValid()) {
00244 edm::LogVerbatim("CkfDebugger") << "RecHit at " << i->recHit()->globalPosition()
00245 << " layer " << ((i->recHit()->det()->geographicalId().rawId() >>16) & 0xF)
00246 << " subdet " << i->recHit()->det()->geographicalId().subdetId()
00247 << " Chi2 " << i->estimate() ;
00248 }
00249 else if (i->recHit()->det() == 0) {
00250 edm::LogVerbatim("CkfDebugger") << "Invalid RecHit returned with zero Det pointer" ;
00251 }
00252 else if (i->recHit()->det() == det(correctHit)) {
00253 edm::LogVerbatim("CkfDebugger") << "Invalid hit returned in correct Det" ;
00254 }
00255 else {
00256 edm::LogVerbatim("CkfDebugger") << "Invalid hit returned in Det at gpos " << i->recHit()->det()->position()
00257 << " correct Det is at " << det(correctHit)->position() ;
00258 }
00259 }
00260
00261
00262 std::pair<CTTRHp, double> correctRecHit =
00263 analyseRecHitExistance( *correctHit, traj.lastMeasurement().updatedState());
00264 if (correctRecHit.first==0 ) {
00265
00266 if ( fabs(correctRecHit.second-0)<0.01 ) {dump[1]++;}
00267 if ( fabs(correctRecHit.second+1)<0.01 ) {dump[8]++;}
00268 if ( fabs(correctRecHit.second+2)<0.01 ) {dump[9]++;}
00269 if ( fabs(correctRecHit.second+3)<0.01 ) {dump[10]++;}
00270 if ( fabs(correctRecHit.second+4)<0.01 ) {dump[11]++;}
00271 if ( fabs(correctRecHit.second+5)<0.01 ) {dump[12]++;}
00272 if ( fabs(correctRecHit.second+6)<0.01 ) {dump[13]++;}
00273 if ( fabs(correctRecHit.second+7)<0.01 ) {dump[14]++;}
00274 if ( fabs(correctRecHit.second+8)<0.01 ) {dump[15]++;}
00275 }
00276 else {
00277
00278 int result = analyseRecHitNotFound(traj,correctRecHit.first);
00279 if (result == 5){
00280 if (correctRecHit.second>30) {
00281 edm::LogVerbatim("CkfDebugger") << "Outling RecHit at pos=" << correctRecHit.first->globalPosition()
00282 << " from SimHit at pos="<< position(correctHit)
00283 << " det=" << correctHit->detUnitId() << " process=" << correctHit->processType() ;
00284 if (hasDelta(correctHit)){
00285 edm::LogVerbatim("CkfDebugger") << "there are deltas on this det" ;
00286 chi2gt30delta++;
00287 dump5[pair<int,int>((correctRecHit.first->det()->geographicalId().subdetId()-1),(layer(correctRecHit.first->det()))-1)]++;
00288 }else{
00289 edm::LogVerbatim("CkfDebugger") << "no deltas on this det" ;
00290 dump[5]++;
00291 chi2gt30++;
00292 dump3[pair<int,int>((correctRecHit.first->det()->geographicalId().subdetId()-1),(layer(correctRecHit.first->det()))-1)]++;
00293 CTTRHp h1 = traj.measurements()[0].recHit();
00294 CTTRHp h2 = traj.measurements()[1].recHit();
00295 TSOS t = traj.measurements()[1].updatedState();
00296 double chi2 = testSeed(h1,h2,t);
00297 if (chi2==-1) {
00298 edm::LogVerbatim("CkfDebugger") << "there were deltas in the seed" ;
00299 chi2gt30deltaSeed++;
00300 }
00301 else {
00302 hchi2seedProb->Fill(chi2);
00303 edm::LogVerbatim("CkfDebugger") << "no deltas in the seed. What is wrong?" ;
00304
00305 TSOS detState = theForwardPropagator->propagate( traj.lastMeasurement().updatedState(), correctRecHit.first->det()->surface());
00306 TSOS simDetState = theForwardPropagator->propagate( traj.lastMeasurement().updatedState(), det(correctHit)->surface());
00307
00308 if (true){
00309 int subdetId = correctRecHit.first->det()->geographicalId().subdetId();
00310 int layerId = layer(correctRecHit.first->det());
00311
00312
00313 LogTrace("CkfDebugger") << "position(correctHit)=" << position(correctHit) ;
00314 LogTrace("CkfDebugger") << "correctRecHit.first->globalPosition()=" << correctRecHit.first->globalPosition() ;
00315 LogTrace("CkfDebugger") << "detState.globalPosition()=" << detState.globalPosition() ;
00316 LogTrace("CkfDebugger") << "simDetState.globalPosition()=" << simDetState.globalPosition() ;
00317
00318 LogTrace("CkfDebugger") << "correctHit->localPosition()=" << correctHit->localPosition() ;
00319 LogTrace("CkfDebugger") << "correctRecHit.first->localPosition()=" << correctRecHit.first->localPosition() ;
00320 LogTrace("CkfDebugger") << "correctRecHit.first->localPositionError()=" << correctRecHit.first->localPositionError() ;
00321 LogTrace("CkfDebugger") << "detState.localPosition()=" << detState.localPosition() ;
00322 LogTrace("CkfDebugger") << "detState.localError().positionError()=" << detState.localError().positionError() ;
00323 LogTrace("CkfDebugger") << "simDetState.localPosition()=" << simDetState.localPosition() ;
00324 LogTrace("CkfDebugger") << "simDetState.localError().positionError()=" << simDetState.localError().positionError() ;
00325 double pullx_shrh = (correctHit->localPosition().x()-correctRecHit.first->localPosition().x())/
00326 sqrt(correctRecHit.first->localPositionError().xx());
00327 double pully_shrh = 0;
00328 if (correctRecHit.first->localPositionError().yy()!=0)
00329 pully_shrh = (correctHit->localPosition().y()-correctRecHit.first->localPosition().y())/
00330 sqrt(correctRecHit.first->localPositionError().yy());
00331 double pullx_shst = (correctHit->localPosition().x()-simDetState.localPosition().x())/
00332 sqrt(simDetState.localError().positionError().xx());
00333 double pully_shst = (correctHit->localPosition().y()-simDetState.localPosition().y())/
00334 sqrt(simDetState.localError().positionError().yy());
00335
00336 LogTrace("CkfDebugger") << "pullx(sh-rh)=" << pullx_shrh ;
00337 LogTrace("CkfDebugger") << "pully(sh-rh)=" << pully_shrh ;
00338 LogTrace("CkfDebugger") << "pullx(sh-st)=" << pullx_shst ;
00339 LogTrace("CkfDebugger") << "pully(sh-st)=" << pully_shst ;
00340
00341 LogTrace("CkfDebugger") << "pullx(st-rh)=" << (detState.localPosition().x()-correctRecHit.first->localPosition().x())/
00342 sqrt(correctRecHit.first->localPositionError().xx()+detState.localError().positionError().xx()) ;
00343
00344 std::pair<double,double> pulls = computePulls(correctRecHit.first, detState);
00345 if (subdetId>0 &&subdetId<7 && layerId>0 && layerId<10) {
00346 stringstream title;
00347 title.str("");
00348 title << "pullX_" << subdetId << "-" << layerId << "_sh-rh";
00349 hPullX_shrh[title.str()]->Fill( pullx_shrh );
00350 title.str("");
00351 title << "pullY_" << subdetId << "-" << layerId << "_sh-rh";
00352 hPullY_shrh[title.str()]->Fill( pully_shrh );
00353 title.str("");
00354 title << "pullX_" << subdetId << "-" << layerId <<"_sh-st";
00355 hPullX_shst[title.str()]->Fill( pullx_shst );
00356 title.str("");
00357 title << "pullY_" << subdetId << "-" << layerId <<"_sh-st";
00358 hPullY_shst[title.str()]->Fill( pully_shst );
00359 title.str("");
00360 title << "pullX_" << subdetId << "-" << layerId <<"_st-rh";
00361 hPullX_strh[title.str()]->Fill(pulls.first);
00362 title.str("");
00363 title << "pullY_" << subdetId << "-" << layerId <<"_st-rh";
00364 hPullY_strh[title.str()]->Fill(pulls.second);
00365
00366 GlobalPoint shGPos = position(correctHit);
00367 GlobalPoint stGPos = simDetState.globalPosition();
00368 GlobalError stGPosErr = simDetState.cartesianError().position();
00369 double pullGPx = (shGPos.x()-stGPos.x())/sqrt(stGPosErr.cxx());
00370 title.str("");
00371 title << "PullGP_X_" << subdetId << "-" << layerId << "_sh-st";
00372 hPullGP_X_shst[title.str()]->Fill(pullGPx);
00373 title.str("");
00374 title << "PullGP_Y_" << subdetId << "-" << layerId << "_sh-st";
00375 hPullGP_Y_shst[title.str()]->Fill((shGPos.y()-stGPos.y())/sqrt(stGPosErr.cyy()));
00376 title.str("");
00377 title << "PullGP_Z_" << subdetId << "-" << layerId << "_sh-st";
00378 hPullGP_Z_shst[title.str()]->Fill((shGPos.z()-stGPos.z())/sqrt(stGPosErr.czz()));
00379
00380 if (subdetId==3&&layerId==1){
00381 hPullGPXvsGPX_shst->Fill(pullGPx,shGPos.x());
00382 hPullGPXvsGPY_shst->Fill(pullGPx,shGPos.y());
00383 hPullGPXvsGPZ_shst->Fill(pullGPx,shGPos.z());
00384 hPullGPXvsGPr_shst->Fill(pullGPx,shGPos.mag());
00385 hPullGPXvsGPeta_shst->Fill(pullGPx,shGPos.eta());
00386 hPullGPXvsGPphi_shst->Fill(pullGPx,shGPos.phi());
00387 }
00388 if (dynamic_cast<const SiStripMatchedRecHit2D*>(correctRecHit.first->hit())) {
00389 LogTrace("CkfDebugger") << "MONO HIT";
00390 CTTRHp tMonoHit = theTTRHBuilder->build(dynamic_cast<const SiStripMatchedRecHit2D*>(correctRecHit.first->hit())->monoHit());
00391 const PSimHit sMonoHit = *(hitAssociator->associateHit(*tMonoHit->hit()).begin());
00392 TSOS monoState = theForwardPropagator->propagate( traj.lastMeasurement().updatedState(), tMonoHit->det()->surface());
00393 double pullM_shrh = (sMonoHit.localPosition().x()-tMonoHit->localPosition().x())/
00394 sqrt(tMonoHit->localPositionError().xx());
00395 double pullM_shst = (sMonoHit.localPosition().x()-monoState.localPosition().x())/
00396 sqrt(monoState.localError().positionError().xx());
00397 std::pair<double,double> pullsMono = computePulls(tMonoHit, monoState);
00398 title.str("");
00399 title << "pullM_" << subdetId << "-" << layerId << "_sh-rh";
00400 hPullM_shrh[title.str()]->Fill(pullM_shrh);
00401 title.str("");
00402 title << "pullM_" << subdetId << "-" << layerId << "_sh-st";
00403 hPullM_shst[title.str()]->Fill(pullM_shst);
00404 title.str("");
00405 title << "pullM_" << subdetId << "-" << layerId << "_st-rh";
00406 hPullM_strh[title.str()]->Fill(pullsMono.first);
00407
00408 LogTrace("CkfDebugger") << "STEREO HIT";
00409 CTTRHp tStereoHit = theTTRHBuilder->build(dynamic_cast<const SiStripMatchedRecHit2D*>(correctRecHit.first->hit())->stereoHit());
00410 const PSimHit sStereoHit = *(hitAssociator->associateHit(*tStereoHit->hit()).begin());
00411 TSOS stereoState = theForwardPropagator->propagate( traj.lastMeasurement().updatedState(), tStereoHit->det()->surface());
00412 double pullS_shrh = (sStereoHit.localPosition().x()-tStereoHit->localPosition().x())/
00413 sqrt(tStereoHit->localPositionError().xx());
00414 double pullS_shst = (sStereoHit.localPosition().x()-stereoState.localPosition().x())/
00415 sqrt(stereoState.localError().positionError().xx());
00416 std::pair<double,double> pullsStereo = computePulls(tStereoHit, stereoState);
00417 title.str("");
00418 title << "pullS_" << subdetId << "-" << layerId << "_sh-rh";
00419 hPullS_shrh[title.str()]->Fill(pullS_shrh);
00420 title.str("");
00421 title << "pullS_" << subdetId << "-" << layerId << "_sh-st";
00422 hPullS_shst[title.str()]->Fill(pullS_shst);
00423 title.str("");
00424 title << "pullS_" << subdetId << "-" << layerId << "_st-rh";
00425 hPullS_strh[title.str()]->Fill(pullsStereo.first);
00426 }
00427 } else
00428 edm::LogVerbatim("CkfDebugger") << "unexpected result: wrong det or layer id "
00429 << subdetId << " " << layerId << " "
00430 << correctRecHit.first->det()->geographicalId().rawId();
00431 }
00432 }
00433 }
00434 }
00435 else {
00436 edm::LogVerbatim("CkfDebugger") << "unexpected result " << correctRecHit.second ;
00437 dump[6]++;chi2ls30++;
00438 }
00439 }
00440 else dump[result]++;
00441 if (result == 3){
00442 dump2[pair<int,int>((correctRecHit.first->det()->geographicalId().subdetId()-1),(layer(correctRecHit.first->det()))-1)]++;
00443 }
00444 if (result == 4){
00445 dump4[pair<int,int>((correctRecHit.first->det()->geographicalId().subdetId()-1),(layer(correctRecHit.first->det()))-1)]++;
00446 }
00447 if (correctRecHit.second>30) {
00448 dump[7]++;totchi2gt30++;
00449 }
00450 }
00451 return false;
00452 }
00453
00454
00455 bool CkfDebugger::correctTrajectory( const Trajectory& traj,unsigned int& trajId) const
00456 {
00457 LogTrace("CkfDebugger") << "now in correctTrajectory" ;
00458 Trajectory::RecHitContainer hits = traj.recHits();
00459
00460 std::vector<SimHitIdpr> currentTrackId = hitAssociator->associateHitId(*hits.front()->hit());
00461 if (currentTrackId.size() == 0) return false;
00462
00463 for (Trajectory::RecHitContainer::const_iterator rh=hits.begin(); rh!=hits.end(); ++rh) {
00464
00465
00466 if (!(*rh)->hit()->isValid()) {
00467
00468 return false;
00469 }
00470
00471
00472 bool nogoodhit = true;
00473 std::vector<PSimHit> assSimHits = hitAssociator->associateHit(*(*rh)->hit());
00474 for (std::vector<PSimHit>::iterator shit=assSimHits.begin();shit!=assSimHits.end();shit++){
00475 if (goodSimHit(*shit)) nogoodhit=false;
00476 }
00477 if (nogoodhit) return false;
00478
00479
00480 bool test = true;
00481 std::vector<SimHitIdpr> nextTrackId = hitAssociator->associateHitId(*(*rh)->hit());
00482 for (std::vector<SimHitIdpr>::iterator i=currentTrackId.begin();i!=currentTrackId.end();i++){
00483 for (std::vector<SimHitIdpr>::iterator j=nextTrackId.begin();j!=nextTrackId.end();j++){
00484 if (i->first == j->first) test = false;
00485
00486 trajId = j->first;
00487 }
00488 }
00489 if (test) {return false;}
00490
00491
00492 }
00493
00494 return true;
00495 }
00496
00497 int CkfDebugger::assocTrackId(CTTRHp rechit) const
00498 {
00499 LogTrace("CkfDebugger") << "now in assocTrackId" ;
00500
00501 if (!rechit->hit()->isValid()) {
00502 return -1;
00503 }
00504
00505 std::vector<SimHitIdpr> ids = hitAssociator->associateHitId(*rechit->hit());
00506 if (ids.size()!=0) {
00507 return ids[0].first;
00508 }
00509 else {
00510 return -1;
00511 }
00512 }
00513
00514
00515 vector<const PSimHit*> CkfDebugger::nextCorrectHits( const Trajectory& traj, unsigned int& trajId)
00516 {
00517 std::vector<const PSimHit*> result;
00518
00519 LogTrace("CkfDebugger") << "now in nextCorrectHits" ;
00520 TransientTrackingRecHit::ConstRecHitPointer lastRecHit = traj.lastMeasurement().recHit();
00521 TransientTrackingRecHit::RecHitContainer comp = lastRecHit->transientHits();
00522 if (!comp.empty()) {
00523 float maxR = 0;
00524 for (TransientTrackingRecHit::RecHitContainer::const_iterator ch=comp.begin();
00525 ch!=comp.end(); ++ch) {
00526 if ((*ch)->globalPosition().mag() > maxR) lastRecHit = *ch;
00527 maxR = (*ch)->globalPosition().mag();
00528 }
00529 }
00530 edm::LogVerbatim("CkfDebugger") << "CkfDebugger: lastRecHit is at gpos " << lastRecHit->globalPosition()
00531 << " layer " << layer((lastRecHit->det()))
00532 << " subdet " << lastRecHit->det()->geographicalId().subdetId() ;
00533
00534
00535 const std::vector<PSimHit>& pSimHitVec = hitAssociator->associateHit(*lastRecHit->hit());
00536 for (std::vector<PSimHit>::const_iterator shit=pSimHitVec.begin();shit!=pSimHitVec.end();shit++){
00537 const GeomDetUnit* detUnit = theTrackerGeom->idToDetUnit( DetId(shit->detUnitId()));
00538 LogTrace("CkfDebugger") << "from hitAssociator SimHits are at GP=" << detUnit->toGlobal( shit->localPosition())
00539 << " traId=" << shit->trackId() << " particleType " << shit->particleType()
00540 << " pabs=" << shit->pabs() << " detUnitId=" << shit->detUnitId() << " layer " << layer((det(&*shit)))
00541 << " subdet " << det(&*shit)->geographicalId().subdetId() ;
00542 }
00543
00544
00545 const PSimHit * lastPSH = 0;
00546 if (!pSimHitVec.empty()) {
00547 float maxTOF = 0;
00548 for (std::vector<PSimHit>::const_iterator ch=pSimHitVec.begin(); ch!=pSimHitVec.end(); ++ch) {
00549 if ( ( ch->trackId()== trajId) && (ch->timeOfFlight() > maxTOF) && ( goodSimHit(*ch) )) {
00550 lastPSH = &*ch;
00551 maxTOF = lastPSH->timeOfFlight();
00552 }
00553 }
00554 }
00555 else return result;
00556 if (lastPSH == 0) return result;
00557 edm::LogVerbatim("CkfDebugger") << "CkfDebugger: corresponding SimHit is at gpos " << position(&*lastPSH) ;
00558
00559
00560 std::vector<PSimHit*> trackHits = idHitsMap[trajId];
00561 if (fabs((double)(trackHits.back()->detUnitId()-lastPSH->detUnitId()))<1 ) return result;
00562 std::vector<PSimHit*>::iterator currentIt = trackHits.end();
00563 for (std::vector<PSimHit*>::iterator it=trackHits.begin();
00564 it!=trackHits.end();it++){
00565 if (goodSimHit(**it) &&
00566 ( lastPSH->timeOfFlight()<(*it)->timeOfFlight() ) &&
00567
00568 ( (det(lastPSH)->geographicalId().subdetId()!=det(*it)->geographicalId().subdetId()) ||
00569 (layer(det(lastPSH))!=layer(det(*it)) ) )
00570 ){
00571 edm::LogVerbatim("CkfDebugger") << "Next good PSimHit is at gpos " << position(*it) ;
00572 result.push_back(*it);
00573 currentIt = it;
00574 break;
00575 }
00576 }
00577 bool samelayer = true;
00578 if (currentIt!=(trackHits.end()-1) && currentIt!=trackHits.end()) {
00579 for (std::vector<PSimHit*>::iterator nextIt = currentIt; (samelayer && nextIt!=trackHits.end()) ;nextIt++){
00580 if (goodSimHit(**nextIt)){
00581 if ( (det(*nextIt)->geographicalId().subdetId()==det(*currentIt)->geographicalId().subdetId()) &&
00582 (layer(det(*nextIt))==layer(det(*currentIt)) ) ) {
00583 result.push_back(*nextIt);
00584 }
00585 else samelayer = false;
00586 }
00587 }
00588 }
00589
00590 return result;
00591 }
00592
00593 bool CkfDebugger::goodSimHit(const PSimHit& sh) const
00594 {
00595 if (sh.pabs() > 0.9) return true;
00596 else return false;
00597 }
00598
00599
00600 bool CkfDebugger::associated(CTTRHp rechit, const PSimHit& pSimHit) const
00601 {
00602 LogTrace("CkfDebugger") << "now in associated" ;
00603
00604 if (!rechit->isValid()) return false;
00605
00606 const std::vector<PSimHit>& pSimHitVec = hitAssociator->associateHit(*rechit->hit());
00607
00608 for (std::vector<PSimHit>::const_iterator shit=pSimHitVec.begin();shit!=pSimHitVec.end();shit++){
00609
00610
00611
00612
00613
00614 if ( ( fabs((*shit).timeOfFlight()-pSimHit.timeOfFlight())<1e-9 ) &&
00615 ( fabs((*shit).pabs()-pSimHit.pabs())<1e-9 ) ) return true;
00616 }
00617 return false;
00618 }
00619
00620 bool CkfDebugger::correctMeas( const TM& tm, const PSimHit* correctHit) const
00621 {
00622 LogTrace("CkfDebugger") << "now in correctMeas" ;
00623 CTTRHp recHit = tm.recHit();
00624 if (recHit->isValid()) LogTrace("CkfDebugger") << "hit at position:" << recHit->globalPosition() ;
00625 TransientTrackingRecHit::RecHitContainer comp = recHit->transientHits();
00626 if (comp.empty()) {
00627
00628 return associated( recHit, *correctHit);
00629 }
00630 else {
00631 for (TransientTrackingRecHit::RecHitContainer::const_iterator ch=comp.begin();
00632 ch!=comp.end(); ++ch) {
00633 if ( associated( recHit, *correctHit)) {
00634
00635 for (TransientTrackingRecHit::RecHitContainer::const_iterator ch2=comp.begin();
00636 ch2!=comp.end(); ++ch2) {
00637 if (ch2 == ch) continue;
00639
00640 bool test=true;
00641 std::vector<SimHitIdpr> ids = hitAssociator->associateHitId(*(*ch2)->hit());
00642 for (std::vector<SimHitIdpr>::iterator j=ids.begin();j!=ids.end();j++){
00643
00644 if (correctHit->trackId()==j->first) {
00645 test=false;
00646
00647 }
00648 }
00649 if (assocTrackId( *ch2) != ((int)( correctHit->trackId())) ) {LogTrace("CkfDebugger") << "returning false 1" ;}
00650 if (test) {
00651
00652 return false;
00653 }
00654
00655
00656
00657 }
00658 return true;
00659 }
00660 }
00661 return false;
00662 }
00663 }
00664
00665
00666 pair<CTTRHp, double> CkfDebugger::analyseRecHitExistance( const PSimHit& sh, const TSOS& startingState)
00667 {
00668 LogTrace("CkfDebugger") << "now in analyseRecHitExistance" ;
00669
00670 std::pair<CTTRHp, double> result;
00671
00672 const MeasurementDet* simHitDet = theMeasurementTracker->idToDet( DetId( sh.detUnitId()));
00673 TSOS simHitState = TSOSFromSimHitFactory()(sh, *det(&sh), *theMagField);
00674 MeasurementDet::RecHitContainer recHits = simHitDet->recHits( simHitState);
00675
00676
00677 TSOS firstDetState = theForwardPropagator->propagate( startingState, det(&sh)->surface());
00678 if (!firstDetState.isValid()) {
00679 edm::LogVerbatim("CkfDebugger") << "CkfDebugger: propagation failed from state " << startingState << " to first det surface "
00680 << position(&sh) ;
00681 propagation++;
00682 return std::pair<CTTRHp, double>((CTTRHp)(0),-1);
00683 }
00684
00685 bool found = false;
00686 for ( MeasurementDet::RecHitContainer::const_iterator rh = recHits.begin(); rh != recHits.end(); rh++) {
00687 if ( associated( *rh, sh)) {
00688 found = true;
00689 result = std::pair<CTTRHp, double>(*rh,theChi2->estimate( firstDetState, **rh).second);
00690 edm::LogVerbatim("CkfDebugger") << "CkfDebugger: A RecHit associated to the correct Simhit exists at lpos "
00691 << (**rh).localPosition()
00692 << " gpos " << (**rh).globalPosition()
00693 << " layer " << layer((**rh).det())
00694 << " subdet " << (**rh).det()->geographicalId().subdetId()
00695 << " Chi2 " << theChi2->estimate( firstDetState, **rh).second;
00696 }
00697 }
00698 if (!found) {
00699 edm::LogVerbatim("CkfDebugger") << "CkfDebugger: there is no RecHit associated to the correct SimHit." ;
00700 edm::LogVerbatim("CkfDebugger") << " There are " << recHits.size() << " RecHits in the simHit DetUnit" ;
00701 edm::LogVerbatim("CkfDebugger") << "SH GP=" << position(&sh) << " subdet=" << det(&sh)->geographicalId().subdetId()
00702 << " layer=" << layer(det(&sh)) ;
00703 int y=0;
00704 for (MeasurementDet::RecHitContainer::const_iterator rh = recHits.begin(); rh != recHits.end(); rh++)
00705 edm::LogVerbatim("CkfDebugger") << "RH#" << y++ << " GP=" << (**rh).globalPosition() << " subdet=" << (**rh).det()->geographicalId().subdetId()
00706 << " layer=" << layer((**rh).det()) ;
00707 for ( MeasurementDet::RecHitContainer::const_iterator rh = recHits.begin(); rh != recHits.end(); rh++) {
00708 edm::LogVerbatim("CkfDebugger") << "Non-associated RecHit at pos " << (**rh).localPosition() ;
00709 }
00710 }
00711
00712 bool found2 = false;
00713 const PSimHit* sh2;
00714 StripSubdetector subdet( det(&sh)->geographicalId());
00715 if (!subdet.glued()) {
00716 edm::LogVerbatim("CkfDebugger") << "The DetUnit is not part of a GluedDet" ;
00717 if (found) {
00718 if (result.second>30){
00719 LogTrace("CkfDebugger") << "rh->parameters()=" << result.first->parameters() ;
00720 LogTrace("CkfDebugger") << "rh->parametersError()=" << result.first->parametersError() ;
00721 MeasurementExtractor me(firstDetState);
00722 AlgebraicVector r(result.first->parameters() - me.measuredParameters(*result.first));
00723 LogTrace("CkfDebugger") << "me.measuredParameters(**rh)=" << me.measuredParameters(*result.first) ;
00724 LogTrace("CkfDebugger") << "me.measuredError(**rh)=" << me.measuredError(*result.first) ;
00725 AlgebraicSymMatrix R(result.first->parametersError() + me.measuredError(*result.first));
00726 LogTrace("CkfDebugger") << "r=" << r ;
00727 LogTrace("CkfDebugger") << "R=" << R ;
00728 int ierr;
00729 R.invert(ierr);
00730 LogTrace("CkfDebugger") << "R(-1)=" << R ;
00731 LogTrace("CkfDebugger") << "chi2=" << R.similarity(r) ;
00732 }
00733 return result;
00734 }
00735 else {
00736 simple_hit_not_found++;
00737 return std::pair<CTTRHp, double>((CTTRHp)(0),-8);
00738 }
00739 } else {
00740 edm::LogVerbatim("CkfDebugger") << "The DetUnit is part of a GluedDet" ;
00741 DetId partnerDetId = DetId( subdet.partnerDetId());
00742
00743 sh2 = pSimHit( sh.trackId(), partnerDetId);
00744 if (sh2 == 0) {
00745 edm::LogVerbatim("CkfDebugger") << "Partner DetUnit does not have a SimHit from the same track" ;
00746 if (found) {
00747
00748 TrackingRecHitProjector<ProjectedRecHit2D> proj;
00749 DetId gid = gluedId( subdet);
00750 const MeasurementDet* gluedDet = theMeasurementTracker->idToDet( gid);
00751 TSOS gluedTSOS = theForwardPropagator->propagate(startingState, gluedDet->geomDet().surface());
00752 CTTRHp projHit = proj.project( *result.first,gluedDet->geomDet(),gluedTSOS).get();
00753
00754
00755 double chi2 = theChi2->estimate(gluedTSOS, *proj.project( *result.first,gluedDet->geomDet(),gluedTSOS)).second;
00756 return std::pair<CTTRHp, double>(projHit,chi2);
00757 }
00758 }
00759 else {
00760 edm::LogVerbatim("CkfDebugger") << "Partner DetUnit has a good SimHit at gpos " << position(sh2)
00761 << " lpos " << sh2->localPosition() ;
00762
00763
00764 const MeasurementDet* partnerDet = theMeasurementTracker->idToDet( partnerDetId);
00765 if (partnerDet == 0) {
00766 edm::LogVerbatim("CkfDebugger") << "Partner measurementDet not found!!!" ;
00767 partner_det_not_fuond++;
00768 return std::pair<CTTRHp, double>((CTTRHp)(0),-3);
00769 }
00770 TSOS simHitState2 = TSOSFromSimHitFactory()(*sh2, *det(sh2), *theMagField);
00771 MeasurementDet::RecHitContainer recHits2 = partnerDet->recHits( simHitState2);
00772
00773 TSOS secondDetState = theForwardPropagator->propagate( startingState, det(sh2)->surface());
00774 if (!secondDetState.isValid()) {
00775 edm::LogVerbatim("CkfDebugger") << "CkfDebugger: propagation failed from state " << startingState << " to second det surface "
00776 << position(sh2) ;
00777 propagation++;
00778 return std::pair<CTTRHp, double>((CTTRHp)(0),-1);
00779 }
00780
00781 for ( MeasurementDet::RecHitContainer::const_iterator rh = recHits2.begin(); rh != recHits2.end(); rh++) {
00782 if ( associated( *rh, *sh2)) {
00783 found2 = true;
00784 edm::LogVerbatim("CkfDebugger") << "CkfDebugger: A RecHit associated to the correct Simhit exists at lpos "
00785 << (**rh).localPosition()
00786 << " gpos " << (**rh).globalPosition()
00787 << " Chi2 " << theChi2->estimate( secondDetState, **rh).second
00788 ;
00789 }
00790 }
00791 if (!found2) {
00792 edm::LogVerbatim("CkfDebugger") << "CkfDebugger: there is no RecHit associated to the correct SimHit." ;
00793 LogTrace("CkfDebugger") << " There are " << recHits.size() << " RecHits in the simHit DetUnit" ;
00794 for ( MeasurementDet::RecHitContainer::const_iterator rh = recHits.begin(); rh != recHits.end(); rh++) {
00795 LogTrace("CkfDebugger") << "Non-associated RecHit at pos " << (**rh).localPosition() ;
00796 }
00797 }
00798 }
00799 }
00800
00801 MeasurementDet::RecHitContainer gluedHits;
00802 if (found && found2) {
00803
00804 DetId gid = gluedId( subdet);
00805 const MeasurementDet* gluedDet = theMeasurementTracker->idToDet( gid);
00806 if ( gluedDet == 0) {
00807 edm::LogVerbatim("CkfDebugger") << "CkfDebugger ERROR: glued MeasurementDet not found!" ;
00808 glued_det_not_fuond++;
00809 return std::pair<CTTRHp, double>((CTTRHp)(0),-4);
00810 }
00811
00812 TSOS gluedDetState = theForwardPropagator->propagate( startingState, gluedDet->surface());
00813 if (!gluedDetState.isValid()) {
00814 edm::LogVerbatim("CkfDebugger") << "CkfDebugger: propagation failed from state " << startingState << " to det surface "
00815 << gluedDet->position() ;
00816 propagation++;
00817 return std::pair<CTTRHp, double>((CTTRHp)(0),-1);
00818 }
00819
00820 gluedHits = gluedDet->recHits( gluedDetState);
00821 edm::LogVerbatim("CkfDebugger") << "CkfDebugger: the GluedDet returned " << gluedHits.size() << " hits" ;
00822 if (gluedHits.size()==0){
00823 edm::LogVerbatim("CkfDebugger") << "Found and associated mono and stereo recHits but not matched!!!" ;
00824 matched_not_found++;
00825 return std::pair<CTTRHp, double>((CTTRHp)(0),-5);
00826 }
00827 bool found3 = false;
00828 for ( MeasurementDet::RecHitContainer::const_iterator rh = gluedHits.begin(); rh != gluedHits.end(); rh++) {
00829 if ( associated( *rh, sh) && associated( *rh, *sh2)) {
00830 double chi2 = theChi2->estimate(gluedDetState, **rh).second;
00831 edm::LogVerbatim("CkfDebugger") << "Matched hit at lpos " << (**rh).localPosition()
00832 << " gpos " << (**rh).globalPosition()
00833 << " has Chi2 " << chi2
00834 ;
00835 result = std::pair<CTTRHp, double>(&**rh,chi2);
00836 found3 = true;
00837 if (chi2>30){
00838 LogTrace("CkfDebugger") << "rh->parameters()=" << (*rh)->parameters() ;
00839 LogTrace("CkfDebugger") << "rh->parametersError()=" << (*rh)->parametersError() ;
00840 MeasurementExtractor me(gluedDetState);
00841 AlgebraicVector r((*rh)->parameters() - me.measuredParameters(**rh));
00842 LogTrace("CkfDebugger") << "me.measuredParameters(**rh)=" << me.measuredParameters(**rh) ;
00843 LogTrace("CkfDebugger") << "me.measuredError(**rh)=" << me.measuredError(**rh) ;
00844 AlgebraicSymMatrix R((*rh)->parametersError() + me.measuredError(**rh));
00845 LogTrace("CkfDebugger") << "r=" << r ;
00846 LogTrace("CkfDebugger") << "R=" << R ;
00847 int ierr;
00848 R.invert(ierr);
00849 LogTrace("CkfDebugger") << "R(-1)=" << R ;
00850 LogTrace("CkfDebugger") << "chi2=" << R.similarity(r) ;
00851 }
00852 break;
00853 }
00854 }
00855 if (found3) return result;
00856 else {
00857 edm::LogVerbatim("CkfDebugger") << "Found and associated mono and stereo recHits. Matched found but not associated!!!" ;
00858 matched_not_associated++;
00859 return std::pair<CTTRHp, double>((CTTRHp)(0),-6);
00860 }
00861 }
00862 else if ( (found && !found2) || (!found && found2) ) {
00863 edm::LogVerbatim("CkfDebugger") << "Only one component is found" ;
00864 only_one_component++;
00865 return std::pair<CTTRHp, double>((CTTRHp)(0),-7);
00866 }
00867 else {
00868 edm::LogVerbatim("CkfDebugger") << "No component is found" ;
00869 no_component++;
00870 return std::pair<CTTRHp, double>((CTTRHp)(0),-2);
00871 }
00872 other++;
00873 return std::pair<CTTRHp, double>((CTTRHp)(0),0);
00874 }
00875
00876 const PSimHit* CkfDebugger::pSimHit(unsigned int tkId, DetId detId)
00877 {
00878 for (std::vector<PSimHit*>::iterator shi=idHitsMap[tkId].begin(); shi!=idHitsMap[tkId].end(); ++shi) {
00879 if ( (*shi)->detUnitId() == detId.rawId() &&
00880
00881 goodSimHit(**shi) ) {
00882 return (*shi);
00883 }
00884 }
00885 return 0;
00886 }
00887
00888 int CkfDebugger::analyseRecHitNotFound(const Trajectory& traj, CTTRHp correctRecHit)
00889 {
00890 unsigned int correctDetId = correctRecHit->det()->geographicalId().rawId();
00891 int correctLayId = layer(correctRecHit->det());
00892 LogTrace("CkfDebugger") << "correct layer id=" << correctLayId ;
00893
00894 TSOS currentState( traj.lastMeasurement().updatedState() );
00895 std::vector<const DetLayer*> nl = traj.lastLayer()->nextLayers( *currentState.freeState(),traj.direction() );
00896 if (nl.empty()) {
00897 edm::LogVerbatim("CkfDebugger") << "no compatible layers" ;
00898 no_layer++;return 2;
00899 }
00900
00901 TkLayerLess lless;
00902 const DetLayer* detLayer = 0;
00903 bool navLayerAfter = false;
00904 bool test = false;
00905 for (std::vector<const DetLayer*>::iterator il = nl.begin(); il != nl.end(); il++) {
00906 if ( dynamic_cast<const PixelBarrelLayer*>(*il) ){
00907 const PixelBarrelLayer* pbl = dynamic_cast<const PixelBarrelLayer*>(*il);
00908 LogTrace("CkfDebugger") << "pbl->specificSurface().bounds().length()=" << pbl->specificSurface().bounds().length() ;
00909 LogTrace("CkfDebugger") << "pbl->specificSurface().bounds().width()=" << pbl->specificSurface().bounds().width() ;
00910 }
00911 int layId = layer(((*(*il)->basicComponents().begin())));
00912 LogTrace("CkfDebugger") << " subdet=" << (*(*il)->basicComponents().begin())->geographicalId().subdetId() << "layer id=" << layId ;
00913 if (layId==correctLayId) {
00914 test = true;
00915 detLayer = &**il;
00916 break;
00917 }
00918 if ( lless( *il, theGeomSearchTracker->detLayer(correctRecHit->det()->geographicalId()) ))
00919 navLayerAfter = true;
00920 }
00921
00922 if (test) {
00923 edm::LogVerbatim("CkfDebugger") << "correct layer taken into account. layer id: " << correctLayId ;
00924 } else if (navLayerAfter){
00925 edm::LogVerbatim("CkfDebugger")<< "SimHit layer after the layers returned by Navigation.";
00926 edm::LogVerbatim("CkfDebugger")<< "Probably a missing SimHit." ;
00927 edm::LogVerbatim("CkfDebugger")<< "check: " << (correctRecHit->det()->geographicalId().subdetId()) << " " << (layer(correctRecHit->det()));
00928 dump6[pair<int,int>((correctRecHit->det()->geographicalId().subdetId()-1),(layer(correctRecHit->det()))-1)]++;
00929 no_sim_hit++;return 16;
00930 }
00931 else {
00932 edm::LogVerbatim("CkfDebugger") << "correct layer NOT taken into account. correct layer id: " << correctLayId ;
00933 layer_not_found++;
00934 return 3;
00935 }
00936
00937 typedef DetLayer::DetWithState DetWithState;
00938 std::vector<DetWithState> compatDets = detLayer->compatibleDets(currentState,*theForwardPropagator,*theChi2);
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948 bool test2 = false;
00949 for (std::vector<DetWithState>::iterator det=compatDets.begin();det!=compatDets.end();det++){
00950 unsigned int detId = det->first->geographicalId().rawId();
00951
00952
00953
00954
00955 if (detId==gluedId(correctRecHit->det()->geographicalId()).rawId()) {
00956 test2=true;
00957 break;
00958 }
00959 }
00960
00961 if (test2){
00962 edm::LogVerbatim("CkfDebugger") << "correct det taken into account. correctDetId is: " << correctDetId
00963 << ". please check chi2." ;
00964 return 5;
00965 }
00966 else {
00967 edm::LogVerbatim("CkfDebugger") << "correct det NOT taken into account. correctDetId: " << correctDetId ;
00968 det_not_found++;return 4;
00969 }
00970
00971 }
00972
00973 double CkfDebugger::testSeed(CTTRHp recHit1, CTTRHp recHit2, TSOS state){
00974
00975
00976 const std::vector<PSimHit>& pSimHitVec1 = hitAssociator->associateHit(*recHit1->hit());
00977 const std::vector<PSimHit>& pSimHitVec2 = hitAssociator->associateHit(*recHit2->hit());
00978
00979 if ( pSimHitVec1.size()==0 || pSimHitVec2.size()==0 || hasDelta(&(*pSimHitVec1.begin())) || hasDelta(&(*pSimHitVec2.begin())) ) {
00980 edm::LogVerbatim("CkfDebugger") << "Seed has delta or problems" ;
00981 return -1;
00982 }
00983
00984
00985
00986
00987
00988
00989
00990
00991 if (pSimHitVec2.size()!=0) {
00992 const PSimHit& simHit = *pSimHitVec2.begin();
00993
00994 double shlp1 = -1/simHit.momentumAtEntry().mag();
00995 double shlp2 = simHit.momentumAtEntry().x()/simHit.momentumAtEntry().z();
00996 double shlp3 = simHit.momentumAtEntry().y()/simHit.momentumAtEntry().z();
00997 double shlp4 = simHit.localPosition().x();
00998 double shlp5 = simHit.localPosition().y();
00999 AlgebraicVector5 v;
01000 v[0] = shlp1;
01001 v[1] = shlp2;
01002 v[2] = shlp3;
01003 v[3] = shlp4;
01004 v[4] = shlp5;
01005
01006
01007
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031 AlgebraicSymMatrix55 R = state.localError().matrix();
01032 R.Invert();
01033 double chi2 = ROOT::Math::Similarity(v-state.localParameters().vector(), R);
01034 LogTrace("CkfDebugger") << "chi2=" << chi2 ;
01035 return chi2;
01036 }
01037
01038 return 0;
01039
01040 }
01041
01042
01043 CkfDebugger::~CkfDebugger(){
01044 for (int it=0; it!=((int)(dump.size())); it++)
01045 edm::LogVerbatim("CkfDebugger") << "dump " << it << " " << dump[it] ;
01046
01047 edm::LogVerbatim("CkfDebugger") ;
01048 edm::LogVerbatim("CkfDebugger") << "seedWithDelta=" << ((double)seedWithDelta/totSeeds) ;
01049 edm::LogVerbatim("CkfDebugger") << "problems=" << ((double)problems/totSeeds) ;
01050 edm::LogVerbatim("CkfDebugger") << "no_sim_hit=" << ((double)no_sim_hit/totSeeds) ;
01051 edm::LogVerbatim("CkfDebugger") << "no_layer=" << ((double)no_layer/totSeeds) ;
01052 edm::LogVerbatim("CkfDebugger") << "layer_not_found=" << ((double)layer_not_found/totSeeds) ;
01053 edm::LogVerbatim("CkfDebugger") << "det_not_found=" << ((double)det_not_found/totSeeds) ;
01054 edm::LogVerbatim("CkfDebugger") << "chi2gt30=" << ((double)chi2gt30/totSeeds) ;
01055 edm::LogVerbatim("CkfDebugger") << "chi2gt30deltaSeed=" << ((double)chi2gt30deltaSeed/totSeeds) ;
01056 edm::LogVerbatim("CkfDebugger") << "chi2gt30delta=" << ((double)chi2gt30delta/totSeeds) ;
01057 edm::LogVerbatim("CkfDebugger") << "chi2ls30=" << ((double)chi2ls30/totSeeds) ;
01058 edm::LogVerbatim("CkfDebugger") << "simple_hit_not_found=" << ((double)simple_hit_not_found/totSeeds) ;
01059 edm::LogVerbatim("CkfDebugger") << "no_component=" << ((double)no_component/totSeeds) ;
01060 edm::LogVerbatim("CkfDebugger") << "only_one_component=" << ((double)only_one_component/totSeeds) ;
01061 edm::LogVerbatim("CkfDebugger") << "matched_not_found=" << ((double)matched_not_found/totSeeds) ;
01062 edm::LogVerbatim("CkfDebugger") << "matched_not_associated=" << ((double)matched_not_associated/totSeeds) ;
01063 edm::LogVerbatim("CkfDebugger") << "partner_det_not_fuond=" << ((double)partner_det_not_fuond/totSeeds) ;
01064 edm::LogVerbatim("CkfDebugger") << "glued_det_not_fuond=" << ((double)glued_det_not_fuond/totSeeds) ;
01065 edm::LogVerbatim("CkfDebugger") << "propagation=" << ((double)propagation/totSeeds) ;
01066 edm::LogVerbatim("CkfDebugger") << "other=" << ((double)other/totSeeds) ;
01067 edm::LogVerbatim("CkfDebugger") << "totchi2gt30=" << ((double)totchi2gt30/totSeeds) ;
01068 edm::LogVerbatim("CkfDebugger") << "totSeeds=" << totSeeds ;
01069 edm::LogVerbatim("CkfDebugger") ;
01070
01071 edm::LogVerbatim("CkfDebugger") << "layer navigation problems:" ;
01072 for (int i=0; i!=6; i++)
01073 for (int j=0; j!=9; j++){
01074 if (i==0 && j>2) break;
01075 if (i==1 && j>1) break;
01076 if (i==2 && j>3) break;
01077 if (i==3 && j>2) break;
01078 if (i==4 && j>5) break;
01079 if (i==5 && j>8) break;
01080 edm::LogVerbatim("CkfDebugger") << "det=" << i+1 << " lay=" << j+1 << " " << dump2[pair<int,int>(i,j)] ;
01081 }
01082 edm::LogVerbatim("CkfDebugger") << "\nlayer with hit having chi2>30:" ;
01083 for (int i=0; i!=6; i++)
01084 for (int j=0; j!=9; j++){
01085 if (i==0 && j>2) break;
01086 if (i==1 && j>1) break;
01087 if (i==2 && j>3) break;
01088 if (i==3 && j>2) break;
01089 if (i==4 && j>5) break;
01090 if (i==5 && j>8) break;
01091 edm::LogVerbatim("CkfDebugger") << "det=" << i+1 << " lay=" << j+1 << " " << dump3[pair<int,int>(i,j)] ;
01092 }
01093 edm::LogVerbatim("CkfDebugger") << "\nlayer with hit having chi2>30 for delta rays:" ;
01094 for (int i=0; i!=6; i++)
01095 for (int j=0; j!=9; j++){
01096 if (i==0 && j>2) break;
01097 if (i==1 && j>1) break;
01098 if (i==2 && j>3) break;
01099 if (i==3 && j>2) break;
01100 if (i==4 && j>5) break;
01101 if (i==5 && j>8) break;
01102 edm::LogVerbatim("CkfDebugger") << "det=" << i+1 << " lay=" << j+1 << " " << dump5[pair<int,int>(i,j)] ;
01103 }
01104 edm::LogVerbatim("CkfDebugger") << "\nlayer with det not found:" ;
01105 for (int i=0; i!=6; i++)
01106 for (int j=0; j!=9; j++){
01107 if (i==0 && j>2) break;
01108 if (i==1 && j>1) break;
01109 if (i==2 && j>3) break;
01110 if (i==3 && j>2) break;
01111 if (i==4 && j>5) break;
01112 if (i==5 && j>8) break;
01113 edm::LogVerbatim("CkfDebugger") << "det=" << i+1 << " lay=" << j+1 << " " << dump4[pair<int,int>(i,j)] ;
01114 }
01115 edm::LogVerbatim("CkfDebugger") << "\nlayer with correct RecHit after missing Sim Hit:" ;
01116 for (int i=0; i!=6; i++)
01117 for (int j=0; j!=9; j++){
01118 if (i==0 && j>2) break;
01119 if (i==1 && j>1) break;
01120 if (i==2 && j>3) break;
01121 if (i==3 && j>2) break;
01122 if (i==4 && j>5) break;
01123 if (i==5 && j>8) break;
01124 edm::LogVerbatim("CkfDebugger") << "det=" << i+1 << " lay=" << j+1 << " " << dump6[pair<int,int>(i,j)] ;
01125 }
01126 hchi2seedAll->Write();
01127 hchi2seedProb->Write();
01128 std::stringstream title;
01129 for (int i=0; i!=6; i++)
01130 for (int j=0; j!=9; j++){
01131 if (i==0 && j>2) break;
01132 if (i==1 && j>1) break;
01133 if (i==2 && j>3) break;
01134 if (i==3 && j>2) break;
01135 if (i==4 && j>5) break;
01136 if (i==5 && j>8) break;
01137 title.str("");
01138 title << "pullX_" << i+1 << "-" << j+1 << "_sh-rh";
01139 hPullX_shrh[title.str()]->Write();
01140 title.str("");
01141 title << "pullY_" << i+1 << "-" << j+1 << "_sh-rh";
01142 hPullY_shrh[title.str()]->Write();
01143 title.str("");
01144 title << "pullX_" << i+1 << "-" << j+1 << "_sh-st";
01145 hPullX_shst[title.str()]->Write();
01146 title.str("");
01147 title << "pullY_" << i+1 << "-" << j+1 << "_sh-st";
01148 hPullY_shst[title.str()]->Write();
01149 title.str("");
01150 title << "pullX_" << i+1 << "-" << j+1 << "_st-rh";
01151 hPullX_strh[title.str()]->Write();
01152 title.str("");
01153 title << "pullY_" << i+1 << "-" << j+1 << "_st-rh";
01154 hPullY_strh[title.str()]->Write();
01155 title.str("");
01156 title << "PullGP_X_" << i+1 << "-" << j+1 << "_sh-st";
01157 hPullGP_X_shst[title.str()]->Write();
01158 title.str("");
01159 title << "PullGP_Y_" << i+1 << "-" << j+1 << "_sh-st";
01160 hPullGP_Y_shst[title.str()]->Write();
01161 title.str("");
01162 title << "PullGP_Z_" << i+1 << "-" << j+1 << "_sh-st";
01163 hPullGP_Z_shst[title.str()]->Write();
01164 if ( ((i==2||i==4)&&(j==0||j==1)) || (i==3||i==5) ){
01165 title.str("");
01166 title << "pullM_" << i+1 << "-" << j+1 << "_sh-rh";
01167 hPullM_shrh[title.str()]->Write();
01168 title.str("");
01169 title << "pullS_" << i+1 << "-" << j+1 << "_sh-rh";
01170 hPullS_shrh[title.str()]->Write();
01171 title.str("");
01172 title << "pullM_" << i+1 << "-" << j+1 << "_sh-st";
01173 hPullM_shst[title.str()]->Write();
01174 title.str("");
01175 title << "pullS_" << i+1 << "-" << j+1 << "_sh-st";
01176 hPullS_shst[title.str()]->Write();
01177 title.str("");
01178 title << "pullM_" << i+1 << "-" << j+1 << "_st-rh";
01179 hPullM_strh[title.str()]->Write();
01180 title.str("");
01181 title << "pullS_" << i+1 << "-" << j+1 << "_st-rh";
01182 hPullS_strh[title.str()]->Write();
01183 }
01184 }
01185 hPullGPXvsGPX_shst->Write();
01186 hPullGPXvsGPY_shst->Write();
01187 hPullGPXvsGPZ_shst->Write();
01188 hPullGPXvsGPr_shst->Write();
01189 hPullGPXvsGPeta_shst->Write();
01190 hPullGPXvsGPphi_shst->Write();
01191
01192
01193 file->Close();
01194 }