00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <vector>
00018 #include <iostream>
00019 #include <cmath>
00020
00021 #include "RoadSearchTrackCandidateMakerAlgorithm.h"
00022
00023 #include "TrackingTools/RoadSearchHitAccess/interface/RoadSearchHitSorting.h"
00024
00025 #include "DataFormats/RoadSearchCloud/interface/RoadSearchCloud.h"
00026 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
00027 #include "DataFormats/Common/interface/OwnVector.h"
00028
00029 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
00030 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00031 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00032 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00033
00034 #include "MagneticField/Engine/interface/MagneticField.h"
00035 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00036
00037 #include "DataFormats/Common/interface/Handle.h"
00038 #include "FWCore/Framework/interface/ESHandle.h"
00039 #include "FWCore/Framework/interface/EventSetup.h"
00040 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00041
00042 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
00043 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00044 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00045 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
00046 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
00047
00048 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00049 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
00050 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
00051 #include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h"
00052 #include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
00053 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
00054 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
00055 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
00056 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
00057
00058 #include "TrackingTools/TrajectoryCleaning/interface/TrajectoryCleaner.h"
00059 #include "TrackingTools/TrajectoryCleaning/interface/TrajectoryCleanerBySharedHits.h"
00060 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
00061 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
00062 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
00063 #include "RecoTracker/TkSeedGenerator/interface/FastLine.h"
00064 #include "TrackingTools/TrackFitters/interface/KFTrajectorySmoother.h"
00065 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
00066 #include "TrackingTools/TrajectoryState/interface/BasicSingleTrajectoryState.h"
00067
00068 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
00069 #include "RecoLocalTracker/SiStripRecHitConverter/interface/SiStripRecHitMatcher.h"
00070 #include "Geometry/TrackerGeometryBuilder/interface/GluedGeomDet.h"
00071
00072 RoadSearchTrackCandidateMakerAlgorithm::RoadSearchTrackCandidateMakerAlgorithm(const edm::ParameterSet& conf) : conf_(conf) {
00073
00074 theNumHitCut = (unsigned int)conf_.getParameter<int>("NumHitCut");
00075 theChi2Cut = conf_.getParameter<double>("HitChi2Cut");
00076
00077 theEstimator = new Chi2MeasurementEstimator(theChi2Cut);
00078 theUpdator = new KFUpdator();
00079 theTrajectoryCleaner = new TrajectoryCleanerBySharedHits;
00080
00081 CosmicReco_ = conf_.getParameter<bool>("StraightLineNoBeamSpotCloud");
00082 CosmicTrackMerging_ = conf_.getParameter<bool>("CosmicTrackMerging");
00083 MinChunkLength_ = conf_.getParameter<int>("MinimumChunkLength");
00084 nFoundMin_ = conf_.getParameter<int>("nFoundMin");
00085
00086 initialVertexErrorXY_ = conf_.getParameter<double>("InitialVertexErrorXY");
00087 splitMatchedHits_ = conf_.getParameter<bool>("SplitMatchedHits");
00088 cosmicSeedPt_ = conf_.getParameter<double>("CosmicSeedPt");
00089
00090 measurementTrackerName_ = conf_.getParameter<std::string>("MeasurementTrackerName");
00091
00092 debug_ = false;
00093 debugCosmics_ = false;
00094
00095 maxPropagationDistance = 1000.0;
00096 }
00097
00098 RoadSearchTrackCandidateMakerAlgorithm::~RoadSearchTrackCandidateMakerAlgorithm() {
00099 delete theEstimator;
00100 delete theUpdator;
00101 delete theTrajectoryCleaner;
00102
00103
00104 }
00105
00106 void RoadSearchTrackCandidateMakerAlgorithm::run(const RoadSearchCloudCollection* input,
00107 const edm::Event& e,
00108 const edm::EventSetup& es,
00109 TrackCandidateCollection &output)
00110 {
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121 edm::ESHandle<TransientTrackingRecHitBuilder> theBuilder;
00122 std::string builderName = conf_.getParameter<std::string>("TTRHBuilder");
00123 es.get<TransientRecHitRecord>().get(builderName,theBuilder);
00124 ttrhBuilder = theBuilder.product();
00125
00126 edm::ESHandle<MeasurementTracker> measurementTrackerHandle;
00127 es.get<CkfComponentsRecord>().get(measurementTrackerName_, measurementTrackerHandle);
00128 theMeasurementTracker = measurementTrackerHandle.product();
00129
00130 std::vector<Trajectory> FinalTrajectories;
00131
00132
00133
00134
00135 edm::ESHandle<TrackerGeometry> tracker;
00136 es.get<TrackerDigiGeometryRecord>().get(tracker);
00137 trackerGeom = tracker.product();
00138
00139 edm::ESHandle<MagneticField> magField_;
00140 es.get<IdealMagneticFieldRecord>().get(magField_);
00141 magField = magField_.product();
00142
00143 NoFieldCosmic_ = (CosmicReco_ && (magField->inTesla(GlobalPoint(0,0,0)).mag() < 0.01));
00144
00145 theMeasurementTracker->update(e);
00146
00147
00148 theAloPropagator = new PropagatorWithMaterial(alongMomentum,.1057,&(*magField));
00149 theRevPropagator = new PropagatorWithMaterial(oppositeToMomentum,.1057,&(*magField));
00150 theAnalyticalPropagator = new AnalyticalPropagator(magField,anyDirection);
00151
00152 thePropagator = theAloPropagator;
00153
00154
00155 theSmoother = new KFTrajectorySmoother(*theRevPropagator, *theUpdator, *theEstimator);
00156
00157
00158 theHitMatcher = new SiStripRecHitMatcher(3.0);
00159
00160
00161
00162
00163 LogDebug("RoadSearch") << "Clean Clouds input size: " << input->size();
00164 if (debug_) std::cout << std::endl << std::endl
00165 << "*** NEW EVENT: Clean Clouds input size: " << input->size() << std::endl;
00166
00167 int i_c = 0;
00168 for ( RoadSearchCloudCollection::const_iterator cloud = input->begin(); cloud != input->end(); ++cloud ) {
00169
00170
00171 RoadSearchCloud::RecHitVector recHits = cloud->recHits();
00172
00173 std::vector<Trajectory> CloudTrajectories;
00174
00175 if (!CosmicReco_){
00176 std::sort(recHits.begin(),recHits.end(),SortHitPointersByGlobalPosition(tracker.product(),alongMomentum));
00177 }
00178 else {
00179 std::sort(recHits.begin(),recHits.end(),SortHitPointersByY(*tracker));
00180 }
00181
00182 const unsigned int nlost_max = 2;
00183
00184
00185
00186 const unsigned int max_layers = 128;
00187
00188
00189 std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > > RecHitsByLayer;
00190 std::map<const DetLayer*, int> cloud_layer_reference;
00191 std::map<const DetLayer*, int>::iterator hiter;
00192 for(RoadSearchCloud::RecHitVector::const_iterator ihit = recHits.begin();
00193 ihit != recHits.end(); ihit++) {
00194
00195 const DetLayer* thisLayer =
00196 theMeasurementTracker->geometricSearchTracker()->detLayer((*ihit)->geographicalId());
00197
00198 std::map<const DetLayer*, int>::const_iterator ilyr = cloud_layer_reference.find(thisLayer);
00199 if (ilyr==cloud_layer_reference.end())
00200 cloud_layer_reference.insert(std::make_pair( thisLayer, RecHitsByLayer.size()));
00201
00202 if (!RecHitsByLayer.empty() && RecHitsByLayer.back().first == thisLayer) {
00203 RecHitsByLayer.back().second.push_back(*ihit);
00204 }
00205 else {
00206 if (ilyr != cloud_layer_reference.end()){
00207 int ilayer = ilyr->second;
00208 (RecHitsByLayer.begin()+ilayer)->second.push_back(*ihit);
00209 }
00210 else{
00211 if (RecHitsByLayer.size() >= max_layers) break;
00212 lstereo[RecHitsByLayer.size()] = false;
00213 if ((*ihit)->localPositionError().yy()<1.) lstereo[RecHitsByLayer.size()] = true;
00214 RoadSearchCloud::RecHitVector rhc;
00215 rhc.push_back(*ihit);
00216 RecHitsByLayer.push_back(std::make_pair(thisLayer, rhc));
00217 }
00218 }
00219 }
00220
00221 LogDebug("RoadSearch")<<"Cloud #"<<i_c<<" has "<<recHits.size()<<" hits in "<<RecHitsByLayer.size()<<" layers ";
00222 if (debug_) std::cout <<"Cloud "<<i_c<<" has "<<recHits.size()<<" hits in " <<RecHitsByLayer.size() << " layers " <<std::endl;;
00223 ++i_c;
00224
00225 if (debug_){
00226 int ntothit = 0;
00227
00228 for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilhv = RecHitsByLayer.begin();
00229 ilhv != RecHitsByLayer.end(); ++ilhv) {
00230 std::cout<<" Layer " << ilhv-RecHitsByLayer.begin() << " has " << ilhv->second.size() << " hits " << std::endl;
00231 }
00232 std::cout<<std::endl;
00233 for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilhv = RecHitsByLayer.begin();
00234 ilhv != RecHitsByLayer.end(); ++ilhv) {
00235 RoadSearchCloud::RecHitVector theLayerHits = ilhv->second;
00236 for (RoadSearchCloud::RecHitVector::const_iterator ihit = theLayerHits.begin();
00237 ihit != theLayerHits.end(); ++ihit) {
00238
00239 GlobalPoint gp = trackerGeom->idToDet((*ihit)->geographicalId())->surface().toGlobal((*ihit)->localPosition());
00240 if (CosmicReco_){
00241 std::cout << " Hit "<< ntothit
00242 << " x/y/z = "
00243 << gp.x() << " " << gp.y() << " " << gp.z()
00244 <<" in layer " << ilhv-RecHitsByLayer.begin()
00245 << " is hit " << (ihit-theLayerHits.begin())+1
00246 << " of " << theLayerHits.size() << std::endl;
00247 }
00248 else {
00249 std::cout << " Hit "<< ntothit
00250 << " r/z = "
00251 << gp.perp() << " " << gp.z()
00252 <<" in layer " << ilhv-RecHitsByLayer.begin()
00253 << " is hit " << (ihit-theLayerHits.begin())+1
00254 << " of " << theLayerHits.size() << std::endl;
00255 }
00256 ntothit++;
00257 }
00258 }
00259 std::cout<<std::endl;
00260 }
00261
00262
00263
00264
00265 for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilyr0 = RecHitsByLayer.begin();
00266 ilyr0 != RecHitsByLayer.end(); ++ilyr0) {
00267
00268 unsigned int ilayer0 = (unsigned int)(ilyr0-RecHitsByLayer.begin());
00269 if (ilayer0 > RecHitsByLayer.size()-MinChunkLength_) continue;
00270
00271 std::vector<Trajectory> ChunkTrajectories;
00272 std::vector<Trajectory> CleanChunks;
00273 bool all_chunk_layers_used = false;
00274
00275 if (debug_) std::cout << "*** START NEW CHUNK --> layer range (" << ilyr0-RecHitsByLayer.begin()
00276 << "-" << RecHitsByLayer.size()-1 << ")";
00277
00278
00279 RoadSearchCloud::RecHitVector recHits_start = ilyr0->second;
00280
00281
00282
00283
00284
00285
00286
00287
00288 std::multimap<int, const DetLayer*> layer_map;
00289 std::map<const DetLayer*, int> layer_reference;
00290
00291 for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilayer = ilyr0+1;
00292 ilayer != RecHitsByLayer.end(); ++ilayer) {
00293 layer_map.insert(std::make_pair(ilayer->second.size(), ilayer->first));
00294 layer_reference.insert(std::make_pair(ilayer->first, ilayer-RecHitsByLayer.begin()));
00295 }
00296
00297 if (debug_) {
00298 std::cout<<std::endl<<" Available layers are: " << std::endl;
00299 for (std::multimap<int, const DetLayer*>::iterator ilm1 = layer_map.begin();
00300 ilm1 != layer_map.end(); ++ilm1) {
00301 std::map<const DetLayer*, int>::iterator ilr = layer_reference.find(ilm1->second);
00302 if (ilr != layer_reference.end() && debug_)
00303 std::cout << "Layer " << ilr->second << " with " << ilm1->first <<" hits" <<std::endl;;
00304 }
00305 }
00306
00307
00308 std::set<const DetLayer*> the_good_layers;
00309 std::vector<const DetLayer*> the_middle_layers;
00310 RoadSearchCloud::RecHitVector the_recHits_middle;
00311
00312
00313 chooseStartingLayers(RecHitsByLayer,ilyr0,layer_map,the_good_layers,the_middle_layers,the_recHits_middle);
00314 if (debug_) {
00315 std::cout << " From new code... With " << the_good_layers.size() << " useful layers: ";
00316 for (std::set<const DetLayer*>::iterator igl = the_good_layers.begin();
00317 igl!= the_good_layers.end(); ++igl){
00318 std::map<const DetLayer*, int>::iterator ilr = layer_reference.find(*igl);
00319 if (ilr != layer_reference.end()) std::cout << " " << ilr->second;
00320 }
00321 std::cout << std::endl;
00322 std::cout << " From new code... and middle layers: ";
00323 for (std::vector<const DetLayer*>::iterator iml = the_middle_layers.begin();
00324 iml!= the_middle_layers.end(); ++iml){
00325 std::map<const DetLayer*, int>::iterator ilr = layer_reference.find(*iml);
00326 if (ilr != layer_reference.end()) std::cout << " " << ilr->second;
00327 }
00328 std::cout << std::endl;
00329 }
00330 RoadSearchCloud::RecHitVector recHits_inner = recHits_start;
00331 RoadSearchCloud::RecHitVector recHits_outer = the_recHits_middle;
00332 std::set<const DetLayer*> good_layers = the_good_layers;
00333 unsigned int ngoodlayers = good_layers.size();
00334
00335 if (debug_)
00336 std::cout<<"Found " << recHits_inner.size() << " inner hits and " << recHits_outer.size() << " outer hits" << std::endl;
00337
00338
00339 std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > > goodHits;
00340
00341 std::set<const DetLayer*> skipped_layers;
00342 std::map<int, const DetLayer*> skipped_layer_detmap;
00343
00344
00345 goodHits.push_back(*ilyr0);
00346
00347 for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilyr = ilyr0+1;
00348 ilyr != RecHitsByLayer.end(); ++ilyr) {
00349 if (good_layers.find(ilyr->first) != good_layers.end()){
00350 goodHits.push_back(*ilyr);
00351 }
00352 else {
00353 skipped_layers.insert(ilyr->first);
00354 std::map<const DetLayer*, int>::iterator ilr = layer_reference.find(ilyr->first);
00355 if (ilr != layer_reference.end())
00356 skipped_layer_detmap.insert(std::make_pair(ilr->second,ilyr->first));
00357 else
00358 if (debug_) std::cout<<"Couldn't find thisLayer to insert into map..."<<std::endl;
00359 }
00360 }
00361
00362
00363 for (RoadSearchCloud::RecHitVector::const_iterator innerHit = recHits_inner.begin();
00364 innerHit != recHits_inner.end(); ++innerHit) {
00365
00366 const DetLayer* innerHitLayer =
00367 theMeasurementTracker->geometricSearchTracker()->detLayer((*innerHit)->geographicalId());
00368
00369
00370 RoadSearchCloud::RecHitVector::iterator middleHit, outerHit;
00371 RoadSearchCloud::RecHitVector::iterator firstHit, lastHit;
00372
00373 bool triplets = (CosmicReco_ && (magField->inTesla(GlobalPoint(0,0,0)).mag() > 0.01));
00374
00375 if (!triplets){
00376 firstHit = recHits_outer.begin();
00377 lastHit = recHits_outer.end();
00378 }
00379 else if (triplets){
00380 firstHit = recHits_outer.begin()+1;
00381 lastHit = recHits_outer.end();
00382 }
00383
00384 for (RoadSearchCloud::RecHitVector::iterator outerHit = firstHit; outerHit != lastHit; ++outerHit) {
00385
00386 const DetLayer* middleHitLayer = 0;
00387 if (triplets){
00388 middleHit = outerHit-1;
00389 middleHitLayer = theMeasurementTracker->geometricSearchTracker()->detLayer((*middleHit)->geographicalId());
00390 }
00391 const DetLayer* outerHitLayer =
00392 theMeasurementTracker->geometricSearchTracker()->detLayer((*outerHit)->geographicalId());
00393 if (middleHitLayer == outerHitLayer) continue;
00394
00395 FreeTrajectoryState fts;
00396 if (!triplets){
00397 if (debug_){
00398 std::map<const DetLayer*, int>::iterator ilro = layer_reference.find(outerHitLayer);
00399 if (ilro != layer_reference.end()) {
00400 std::cout << "Try trajectory with Inner Hit on Layer " << ilayer0 << " and " ;
00401 std::cout << "Outer Hit on Layer " << ilro->second << std::endl;
00402 }
00403 }
00404 fts = initialTrajectory(es,*innerHit,*outerHit);
00405 }
00406 else if (triplets){
00407 if (debug_){
00408 std::map<const DetLayer*, int>::iterator ilrm = layer_reference.find(middleHitLayer);
00409 std::map<const DetLayer*, int>::iterator ilro = layer_reference.find(outerHitLayer);
00410 if (ilro != layer_reference.end() && ilrm != layer_reference.end()) {
00411 std::cout << "Try trajectory with Hits on Layers " << ilayer0 << " , "
00412 << ilrm->second << " and " << ilro->second << std::endl;
00413 }
00414 }
00415 fts = initialTrajectoryFromTriplet(es,*innerHit,*middleHit,*outerHit);
00416 }
00417
00418 if (!fts.hasError()) continue;
00419 if (debug_) std::cout<<"FTS: " << fts << std::endl;
00420
00421 Trajectory seedTraj = createSeedTrajectory(fts,*innerHit,innerHitLayer);
00422
00423 std::vector<Trajectory> rawTrajectories;
00424 if (seedTraj.isValid() && !seedTraj.measurements().empty() ) rawTrajectories.push_back(seedTraj);
00425
00426
00427 int ntested = 0;
00428
00429 std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilyr_start = (goodHits.begin()+1);
00430 for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilhv = ilyr_start;
00431 ilhv != goodHits.end(); ++ilhv) {
00432 RoadSearchCloud::RecHitVector& hits = ilhv->second;
00433
00434 ++ntested;
00435 if (debug_){
00436 std::map<const DetLayer*, int>::iterator ilr = cloud_layer_reference.find(ilhv->first);
00437 if (ilr != cloud_layer_reference.end())
00438 std::cout << "extrapolating " << rawTrajectories.size()
00439 << " trajectories to layer " << ilr->second
00440 << " which has " << hits.size() << " hits " << std::endl;
00441 }
00442
00443 std::vector<Trajectory>newTrajectories;
00444 for (std::vector<Trajectory>::const_iterator it = rawTrajectories.begin();
00445 it != rawTrajectories.end(); it++) {
00446 if (debug_) std::cout << "extrapolating Trajectory #" << it-rawTrajectories.begin() << std::endl;
00447 if (it->direction()==alongMomentum) thePropagator = theAloPropagator;
00448 else thePropagator = theRevPropagator;
00449
00450 std::vector<Trajectory> theTrajectories = extrapolateTrajectory(*it,hits,
00451 innerHitLayer, *outerHit, outerHitLayer);
00452 if (theTrajectories.empty()) {
00453 if (debug_) std::cout<<" Could not add the hit in this layer " << std::endl;
00454 if (debug_){
00455 std::cout << " --> trajectory " << it-rawTrajectories.begin()
00456 << " has "<<it->recHits().size()<<" hits after "
00457 << (ilhv-ilyr_start+1) << " tested (ntested=" <<ntested<<") "
00458 << " --> misses="<< (ilhv-ilyr_start+1)-(it->recHits().size()-1)
00459 << " but there are " << (goodHits.end() - ilhv)
00460 <<" more layers in first pass and "<< skipped_layers.size() <<" skipped layers " <<std::endl;
00461
00462
00463 }
00464
00465 if ((ilhv-ilyr_start+1)-(it->recHits().size()-1) <= nlost_max){
00466 newTrajectories.push_back(*it);
00467 }
00468 }
00469 else{
00470 for (std::vector<Trajectory>::const_iterator it = theTrajectories.begin();
00471 it != theTrajectories.end(); it++) {
00472 newTrajectories.push_back(*it);
00473 }
00474 }
00475 }
00476 rawTrajectories = newTrajectories;
00477 if (newTrajectories.empty()) break;
00478 }
00479 if (rawTrajectories.size()==0){
00480 continue;
00481 if (debug_) std::cout<<" --> yields ZERO raw trajectories!" << std::endl;
00482 }
00483 if (debug_){
00484 for (std::vector<Trajectory>::const_iterator it = rawTrajectories.begin();
00485 it != rawTrajectories.end(); it++) {
00486 std::cout << " --> yields trajectory with "<<it->recHits().size()<<" hits with chi2="
00487 <<it->chiSquared()<<" and is valid? "<<it->isValid() <<std::endl;
00488 }
00489 }
00490 std::vector<Trajectory> rawCleaned;
00491 theTrajectoryCleaner->clean(rawTrajectories);
00492 for (std::vector<Trajectory>::const_iterator itr = rawTrajectories.begin();
00493 itr != rawTrajectories.end(); ++itr) {
00494
00495 if (!itr->isValid()) continue;
00496 std::set<const DetLayer*> used_layers;
00497 Trajectory::DataContainer tmv = itr->measurements();
00498 for (Trajectory::DataContainer::iterator itm = tmv.begin();
00499 itm != tmv.end(); ++itm) {
00500 TransientTrackingRecHit::ConstRecHitPointer rh = itm->recHit();
00501 if (!rh->isValid()) continue;
00502 used_layers.insert(theMeasurementTracker->geometricSearchTracker()->detLayer(rh->geographicalId()));
00503 }
00504
00505
00506 if (debug_) std::cout<<"Used " << (used_layers.size()-1) << " layers out of " << ngoodlayers
00507 << " good layers, so " << ngoodlayers - (used_layers.size()-1) << " missed "
00508 << std::endl;
00509 if ((int)used_layers.size() < nFoundMin_) continue;
00510 unsigned int nlostlayers = ngoodlayers - (used_layers.size()-1);
00511 if (nlostlayers > nlost_max) continue;
00512
00513 rawCleaned.push_back(*itr);
00514
00515 }
00516 if (!rawCleaned.empty()) {
00517 ChunkTrajectories.insert(ChunkTrajectories.end(), rawCleaned.begin(), rawCleaned.end());
00518 }
00519 }
00520 }
00521
00522
00523
00524
00525 if (debug_) std::cout << "Clean the " << ChunkTrajectories.size()<<" trajectories for this chunk" << std::endl;
00526
00527 theTrajectoryCleaner->clean(ChunkTrajectories);
00528 for (std::vector<Trajectory>::const_iterator it = ChunkTrajectories.begin();
00529 it != ChunkTrajectories.end(); it++) {
00530 if (it->isValid()) CleanChunks.push_back(*it);
00531 }
00532 if (debug_) std::cout <<"After cleaning there are " << CleanChunks.size() << " trajectories for this chunk" << std::endl;
00533
00534
00535
00536
00537
00538
00539
00540
00541 std::vector<Trajectory> extendedChunks;
00542
00543
00544
00545
00546 if (debug_){
00547 if (skipped_layers.empty()) {
00548 std::cout << "all layers were used in first pass" << std::endl;
00549 } else {
00550 std::cout << "There are " << skipped_layer_detmap.size() << " skipped layers:";
00551 for (std::map<int, const DetLayer*>::const_iterator imap = skipped_layer_detmap.begin();
00552 imap!=skipped_layer_detmap.end(); imap++){
00553 std::cout<< " " <<imap->first;
00554 }
00555 std::cout << std::endl;
00556 }
00557 }
00558
00559 for (std::vector<Trajectory>::const_iterator i = CleanChunks.begin();
00560 i != CleanChunks.end(); i++) {
00561 if (!(*i).isValid()) continue;
00562 if (debug_) std::cout<< "Now process CleanChunk trajectory " << i-CleanChunks.begin() << std::endl;
00563 bool all_layers_used = false;
00564 if (skipped_layers.empty() && i->measurements().size() >= theNumHitCut) {
00565 if (debug_) std::cout<<"The trajectory has " << i->measurements().size()
00566 << " hits from a cloud of " << RecHitsByLayer.size()
00567 << " layers and a chunk of " << (RecHitsByLayer.size() - ilayer0) << " layers " << std::endl;
00568 extendedChunks.insert(extendedChunks.end(), *i);
00569 if (i->measurements().size() >= (RecHitsByLayer.size() - ilayer0)){
00570 all_layers_used = true;
00571 break;
00572 }
00573 }
00574 else {
00575
00576 Trajectory temWptraj = *i;
00577 Trajectory::DataContainer tmv = (*i).measurements();
00578 if (tmv.size()+skipped_layer_detmap.size() < theNumHitCut) continue;
00579
00580
00581 if (debug_){
00582 for (Trajectory::DataContainer::const_iterator ih=tmv.begin();
00583 ih!=tmv.end();++ih){
00584 TransientTrackingRecHit::ConstRecHitPointer rh = ih->recHit();
00585 const DetLayer* Layer =
00586 theMeasurementTracker->geometricSearchTracker()->detLayer(rh->geographicalId());
00587 std::map<const DetLayer*, int>::iterator ilr = cloud_layer_reference.find(Layer);
00588 if (ilr != cloud_layer_reference.end())
00589 std::cout << " Hit #"<<ih-tmv.begin() << " of " << tmv.size()
00590 <<" is on Layer " << ilr->second << std::endl;
00591 else
00592 std::cout << " Layer for Hit #"<<ih-tmv.begin() <<" can't be found " << std::endl;
00593 std::cout<<" updatedState:\n" << ih->updatedState() << std::endl;
00594 std::cout<<" predictState:\n" << ih->predictedState() << std::endl;
00595 }
00596 }
00597
00598
00599
00600 std::set<const DetLayer*> final_layers;
00601 Trajectory::DataContainer::const_iterator im = tmv.begin();
00602 Trajectory::DataContainer::const_iterator im2 = tmv.begin();
00603
00604 TrajectoryMeasurement firstMeasurement = i->firstMeasurement();
00605 const DetLayer* firstDetLayer =
00606 theMeasurementTracker->geometricSearchTracker()->detLayer(firstMeasurement.recHit()->geographicalId());
00607
00608 std::vector<Trajectory> freshStartv = theSmoother->trajectories(*i);
00609
00610 if(debug_) std::cout<<"Smoothing has returned " << freshStartv.size() <<" trajectory " << std::endl;
00611 if (!freshStartv.empty()){
00612 if(debug_) std::cout<< "Smoothing of trajectory " <<i-CleanChunks.begin() << " has succeeded with " << freshStartv.begin()->measurements().size() << " hits and a chi2 of " << freshStartv.begin()->chiSquared() <<" for " << freshStartv.begin()->ndof() << " DOF. Now add hits." <<std::endl;
00613 } else {
00614 if (debug_) std::cout<< "Smoothing of trajectory " <<i-CleanChunks.begin() <<" has failed"<<std::endl;
00615 continue;
00616 }
00617
00618 Trajectory freshStart = *freshStartv.begin();
00619 std::vector<TrajectoryMeasurement> freshStartTM = freshStart.measurements();
00620
00621 if (debug_) {
00622 for (std::vector<TrajectoryMeasurement>::const_iterator itm = freshStartTM.begin();itm != freshStartTM.end(); ++itm){
00623 std::cout<<"Trajectory hit " << itm-freshStartTM.begin() <<" updatedState:\n" << itm->updatedState() << std::endl;
00624 }
00625 }
00626
00627 TrajectoryStateOnSurface NewFirstTsos = freshStart.lastMeasurement().updatedState();
00628 if(debug_) std::cout<<"NewFirstTSOS is valid? " << NewFirstTsos.isValid() << std::endl;
00629 if(debug_) std::cout << " NewFirstTSOS:\n " << NewFirstTsos << std::endl;
00630 TransientTrackingRecHit::ConstRecHitPointer rh = freshStart.lastMeasurement().recHit();
00631
00632 if(debug_) {
00633 std::cout<< "First hit for fresh start on det " << rh->det() << ", r/phi/z = " << rh->globalPosition().perp() << " " << rh->globalPosition().phi() << " " << rh->globalPosition().z();
00634 }
00635
00636 PTrajectoryStateOnDet const & pFirstState = trajectoryStateTransform::persistentState(NewFirstTsos,
00637 rh->geographicalId().rawId());
00638 edm::OwnVector<TrackingRecHit> newHits;
00639 newHits.push_back(rh->hit()->clone());
00640
00641 TrajectorySeed tmpseed = TrajectorySeed(pFirstState,
00642 newHits,
00643 i->direction());
00644
00645 thePropagator = theAloPropagator;
00646 if (i->direction()==oppositeToMomentum) thePropagator = theRevPropagator;
00647
00648
00649
00650 Trajectory newTrajectory(tmpseed,tmpseed.direction());
00651
00652 const GeomDet* det = trackerGeom->idToDet(rh->geographicalId());
00653 TrajectoryStateOnSurface invalidState(new BasicSingleTrajectoryState(det->surface()));
00654 newTrajectory.push(TrajectoryMeasurement(invalidState, NewFirstTsos, rh, 0, firstDetLayer));
00655 final_layers.insert(firstDetLayer);
00656
00657 if(debug_) std::cout << "TRAJ is valid: " << newTrajectory.isValid() <<std::endl;
00658
00659 TrajectoryStateOnSurface testTsos = newTrajectory.measurements().back().updatedState();
00660
00661 if(debug_) {
00662 std::cout << "testTSOS is valid!!" << testTsos.isValid() << std::endl;
00663 std::cout << " testTSOS (x/y/z): " << testTsos.globalPosition().x()<< " / " << testTsos.globalPosition().y()<< " / " << testTsos.globalPosition().z() << std::endl;
00664 std::cout << " local position x: " << testTsos.localPosition().x() << "+-" << sqrt(testTsos.localError().positionError().xx()) << std::endl;
00665 }
00666
00667 if (firstDetLayer != ilyr0->first){
00668 if (debug_) std::cout<<"!!! ERROR !!! firstDetLayer ("<<firstDetLayer<<") != ilyr0 ( " <<ilyr0->first <<")"<< std::endl;
00669 }
00670 ++im;
00671
00672 if (debug_){
00673 std::map<const DetLayer*, int>::iterator ilr = cloud_layer_reference.find(firstDetLayer);
00674 if (ilr != cloud_layer_reference.end() ){
00675 std::cout << " First hit is on layer " << ilr->second << std::endl;
00676 }
00677 }
00678 for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilyr = ilyr0+1;
00679 ilyr != RecHitsByLayer.end(); ++ilyr) {
00680
00681 TrajectoryStateOnSurface predTsos;
00682 TrajectoryStateOnSurface currTsos;
00683 TrajectoryMeasurement tm;
00684
00685 if(debug_)
00686 std::cout<<"Trajectory has " << newTrajectory.measurements().size() << " measurements with " << (RecHitsByLayer.end()-ilyr)
00687 << " remaining layers " << std::endl;
00688
00689 if (im != tmv.end()) im2 = im;
00690 TransientTrackingRecHit::ConstRecHitPointer rh = im2->recHit();
00691 if (rh->isValid() &&
00692 theMeasurementTracker->geometricSearchTracker()->detLayer(rh->geographicalId()) == ilyr->first) {
00693
00694 if (debug_) std::cout<<" Layer " << ilyr-RecHitsByLayer.begin() <<" has a good hit " << std::endl;
00695 ++im;
00696
00697 const GeomDet* det = trackerGeom->idToDet(rh->geographicalId());
00698
00699 currTsos = newTrajectory.measurements().back().updatedState();
00700 predTsos = thePropagator->propagate(currTsos, det->surface());
00701 if (!predTsos.isValid()) continue;
00702 GlobalVector propagationDistance = predTsos.globalPosition() - currTsos.globalPosition();
00703 if (propagationDistance.mag() > maxPropagationDistance) continue;
00704 MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rh);
00705 if(debug_) {
00706 std::cout << "Propagation distance2 is " << propagationDistance.mag() << std::endl;
00707 std::cout << "predTSOS is valid!!" << std::endl;
00708 std::cout << " predTSOS (x/y/z): " << predTsos.globalPosition().x()<< " / " << predTsos.globalPosition().y()<< " / " << predTsos.globalPosition().z() << std::endl;
00709 std::cout << " local position x: " << predTsos.localPosition().x() << "+-" << sqrt(predTsos.localError().positionError().xx()) << std::endl;
00710 std::cout << " local position y: " << predTsos.localPosition().y() << "+-" << sqrt(predTsos.localError().positionError().yy()) << std::endl;
00711 std::cout << "currTSOS is valid!! " << currTsos.isValid() << std::endl;
00712 std::cout << " currTSOS (x/y/z): " << currTsos.globalPosition().x()<< " / " << currTsos.globalPosition().y()<< " / " << currTsos.globalPosition().z() << std::endl;
00713 std::cout << " local position x: " << currTsos.localPosition().x() << "+-" << sqrt(currTsos.localError().positionError().xx()) << std::endl;
00714 std::cout << " local position y: " << currTsos.localPosition().y() << "+-" << sqrt(currTsos.localError().positionError().yy()) << std::endl;
00715 }
00716
00717 if (!est.first) {
00718 if (debug_) std::cout<<"Failed to add one of the original hits on a low occupancy layer!!!!" << std::endl;
00719 continue;
00720 }
00721 currTsos = theUpdator->update(predTsos, *rh);
00722 tm = TrajectoryMeasurement(predTsos, currTsos, &(*rh),est.second,ilyr->first);
00723
00724 const TrajectoryStateOnSurface theTSOS = newTrajectory.lastMeasurement().updatedState();
00725
00726
00727
00728
00729
00730
00731
00732 newTrajectory.push(tm,est.second);
00733 final_layers.insert(ilyr->first);
00734 }
00735 else{
00736 if (debug_) std::cout<<" Layer " << ilyr-RecHitsByLayer.begin() <<" is one of the skipped layers " << std::endl;
00737
00738
00739 edm::OwnVector<TrackingRecHit> skipped_hits;
00740 std::set<const GeomDet*> dets;
00741 for (RoadSearchCloud::RecHitVector::const_iterator ih = ilyr->second.begin();
00742 ih != ilyr->second.end(); ++ih) {
00743 skipped_hits.push_back((*ih)->clone());
00744 dets.insert(trackerGeom->idToDet((*ih)->geographicalId()));
00745 }
00746
00747 if(debug_){
00748 std::cout<<" ---> probing missing hits (nh="<< skipped_hits.size() << ", nd=" << dets.size()
00749 << ") in layer " << ilyr-RecHitsByLayer.begin() <<std::endl;
00750 }
00751
00752 const TrajectoryStateOnSurface theTSOS = newTrajectory.lastMeasurement().updatedState();
00753
00754
00755
00756
00757
00758 std::vector<TrajectoryMeasurement> theGoodHits = FindBestHits(theTSOS,dets,theHitMatcher,skipped_hits);
00759 if (!theGoodHits.empty()){
00760 final_layers.insert(ilyr->first);
00761 if (debug_) std::cout<<"Found " << theGoodHits.size() << " good hits to add" << std::endl;
00762 for (std::vector<TrajectoryMeasurement>::const_iterator im=theGoodHits.begin();im!=theGoodHits.end();++im){
00763 newTrajectory.push(*im,im->estimate());
00764 }
00765 }
00766 }
00767 }
00768
00769 if (debug_) std::cout<<"Finished loop over layers in cloud. Trajectory now has " <<newTrajectory.measurements().size()
00770 << " hits. " << std::endl;
00771 if (debug_) std::cout<<"The trajectory has " << newTrajectory.measurements().size() <<" hits on " << final_layers.size()
00772 << " layers from a cloud of " << RecHitsByLayer.size()
00773 << " layers and a chunk of " << (RecHitsByLayer.size() - ilayer0) << " layers " << std::endl;
00774 if (newTrajectory.measurements().size() >= theNumHitCut)
00775 extendedChunks.insert(extendedChunks.end(), newTrajectory);
00776 if (final_layers.size() >= (RecHitsByLayer.size() - ilayer0)){
00777 if (debug_) std::cout<<"All layers of the chunk have been used..." << std::endl;
00778 all_layers_used = true;
00779 }
00780 }
00781 if (all_layers_used) {
00782 if (debug_) std::cout << "All layers were used, so break....." << std::endl;
00783 all_chunk_layers_used = true;
00784 break;
00785 }
00786 if (debug_) std::cout<<"Going to next clean chunk....." << std::endl;
00787 }
00788 if (debug_) std::cout<< "Now Clean the " << extendedChunks.size() << " extended chunks " <<std::endl;
00789 if (extendedChunks.size() > 1) theTrajectoryCleaner->clean(extendedChunks);
00790 for (std::vector<Trajectory>::const_iterator it = extendedChunks.begin();
00791 it != extendedChunks.end(); it++) {
00792 if (it->isValid()) CloudTrajectories.push_back(*it);
00793 }
00794 if (all_chunk_layers_used) break;
00795 }
00796
00797
00798
00799 if (debug_) std::cout<< "Finished with the cloud, so clean the "
00800 << CloudTrajectories.size() << " cloud trajectories "<<std::endl ;
00801 theTrajectoryCleaner->clean(CloudTrajectories);
00802 for (std::vector<Trajectory>::const_iterator it = CloudTrajectories.begin();
00803 it != CloudTrajectories.end(); it++) {
00804 if (it->isValid()) FinalTrajectories.push_back(*it);
00805 }
00806
00807 if (debug_) std::cout<<" Final trajectories now has size " << FinalTrajectories.size()<<std::endl ;
00808
00809 }
00810
00811
00812 if (debug_) std::cout<< " Finished loop over all clouds " <<std::endl;
00813
00814 output = PrepareTrackCandidates(FinalTrajectories);
00815
00816 delete theAloPropagator;
00817 delete theRevPropagator;
00818 delete theAnalyticalPropagator;
00819 delete theHitMatcher;
00820 delete theSmoother;
00821
00822 if (debug_ || debugCosmics_) std::cout<< "Found " << output.size() << " track candidate(s)."<<std::endl;
00823
00824 }
00825
00826
00827
00828 std::vector<TrajectoryMeasurement>
00829 RoadSearchTrackCandidateMakerAlgorithm::FindBestHitsByDet(const TrajectoryStateOnSurface& tsosBefore,
00830 const std::set<const GeomDet*>& theDets,
00831 edm::OwnVector<TrackingRecHit>& theHits)
00832
00833 {
00834
00835
00836 std::vector<TrajectoryMeasurement> theBestHits;
00837
00838
00839 std::map<const GeomDet*, TrajectoryStateOnSurface> dmmap;
00840 for (std::set<const GeomDet*>::iterator idet = theDets.begin();
00841 idet != theDets.end(); ++idet) {
00842 TrajectoryStateOnSurface predTsos = thePropagator->propagate(tsosBefore, (**idet).surface());
00843 if (predTsos.isValid()) {
00844 GlobalVector propagationDistance = predTsos.globalPosition() - tsosBefore.globalPosition();
00845 if (propagationDistance.mag() > maxPropagationDistance) continue;
00846 dmmap.insert(std::make_pair(*idet, predTsos));
00847 }
00848 }
00849
00850 std::map<const GeomDet*, TrajectoryMeasurement> dtmmap;
00851 for (edm::OwnVector<TrackingRecHit>::const_iterator ih = theHits.begin();
00852 ih != theHits.end(); ++ih) {
00853 const GeomDet* det = trackerGeom->idToDet(ih->geographicalId());
00854
00855 std::map<const GeomDet*, TrajectoryStateOnSurface>::iterator idm = dmmap.find(det);
00856 if (idm == dmmap.end()) continue;
00857 TrajectoryStateOnSurface predTsos = idm->second;
00858 TransientTrackingRecHit::RecHitPointer rhit = ttrhBuilder->build(&(*ih));
00859 MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rhit);
00860
00861
00862 if (est.first) {
00863 TrajectoryStateOnSurface currTsos = theUpdator->update(predTsos, *rhit);
00864 std::map<const GeomDet*, TrajectoryMeasurement>::iterator idtm = dtmmap.find(det);
00865 if (idtm == dtmmap.end()) {
00866 TrajectoryMeasurement tm(predTsos, currTsos, &(*rhit),est.second,
00867 theMeasurementTracker->geometricSearchTracker()->detLayer(ih->geographicalId()));
00868 dtmmap.insert(std::make_pair(det, tm));
00869 } else if (idtm->second.estimate() > est.second) {
00870 dtmmap.erase(idtm);
00871 TrajectoryMeasurement tm(predTsos, currTsos, &(*rhit),est.second,
00872 theMeasurementTracker->geometricSearchTracker()->detLayer(ih->geographicalId()));
00873 dtmmap.insert(std::make_pair(det, tm));
00874 }
00875 }
00876 }
00877
00878 if (!dtmmap.empty()) {
00879 for (std::map<const GeomDet*, TrajectoryMeasurement>::iterator idtm = dtmmap.begin();
00880 idtm != dtmmap.end(); ++idtm) {
00881 TrajectoryMeasurement itm = idtm->second;
00882 theBestHits.push_back(itm);
00883 }
00884 }
00885
00886 return theBestHits;
00887 }
00888
00889
00890 std::vector<TrajectoryMeasurement>
00891 RoadSearchTrackCandidateMakerAlgorithm::FindBestHit(const TrajectoryStateOnSurface& tsosBefore,
00892 const std::set<const GeomDet*>& theDets,
00893 edm::OwnVector<TrackingRecHit>& theHits)
00894 {
00895
00896 std::vector<TrajectoryMeasurement> theBestHits;
00897
00898 double bestchi = 10000.0;
00899
00900 std::map<const GeomDet*, TrajectoryStateOnSurface> dmmap;
00901 for (std::set<const GeomDet*>::iterator idet = theDets.begin();
00902 idet != theDets.end(); ++idet) {
00903 TrajectoryStateOnSurface predTsos = thePropagator->propagate(tsosBefore, (**idet).surface());
00904 if (predTsos.isValid()) {
00905 GlobalVector propagationDistance = predTsos.globalPosition() - tsosBefore.globalPosition();
00906 if (propagationDistance.mag() > maxPropagationDistance) continue;
00907 dmmap.insert(std::make_pair(*idet, predTsos));
00908 }
00909 }
00910
00911 std::map<const GeomDet*, TrajectoryMeasurement> dtmmap;
00912 for (edm::OwnVector<TrackingRecHit>::const_iterator ih = theHits.begin();
00913 ih != theHits.end(); ++ih) {
00914 const GeomDet* det = trackerGeom->idToDet(ih->geographicalId());
00915
00916 std::map<const GeomDet*, TrajectoryStateOnSurface>::iterator idm = dmmap.find(det);
00917 if (idm == dmmap.end()) continue;
00918 TrajectoryStateOnSurface predTsos = idm->second;
00919 TransientTrackingRecHit::RecHitPointer rhit = ttrhBuilder->build(&(*ih));
00920 MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rhit);
00921
00922
00923 if (est.first) {
00924 TrajectoryStateOnSurface currTsos = theUpdator->update(predTsos, *rhit);
00925 if (est.second < bestchi){
00926 if(!theBestHits.empty()){
00927 theBestHits.erase(theBestHits.begin());
00928 }
00929 bestchi = est.second;
00930 TrajectoryMeasurement tm(predTsos, currTsos, &(*rhit),est.second,
00931 theMeasurementTracker->geometricSearchTracker()->detLayer(ih->geographicalId()));
00932 theBestHits.push_back(tm);
00933 }
00934 }
00935 }
00936
00937 if (theBestHits.empty()){
00938 if (debug_) std::cout<< "no hits to add! " <<std::endl;
00939 }
00940 else{
00941 for (std::vector<TrajectoryMeasurement>::const_iterator im=theBestHits.begin();im!=theBestHits.end();++im)
00942 if (debug_) std::cout<<" Measurement on layer "
00943 << theMeasurementTracker->geometricSearchTracker()->detLayer(im->recHit()->geographicalId())
00944 << " with estimate " << im->estimate()<<std::endl ;
00945 }
00946
00947 return theBestHits;
00948 }
00949
00950 std::vector<TrajectoryMeasurement>
00951 RoadSearchTrackCandidateMakerAlgorithm::FindBestHits(const TrajectoryStateOnSurface& tsosBefore,
00952 const std::set<const GeomDet*>& theDets,
00953 const SiStripRecHitMatcher* theHitMatcher,
00954 edm::OwnVector<TrackingRecHit>& theHits)
00955
00956 {
00957
00958
00959 std::vector<TrajectoryMeasurement> theBestHits;
00960
00961 TrajectoryMeasurement theBestTM;
00962 bool firstTM = true;
00963
00964
00965 std::map<const GeomDet*, TrajectoryStateOnSurface> dmmap;
00966 for (std::set<const GeomDet*>::iterator idet = theDets.begin();
00967 idet != theDets.end(); ++idet) {
00968 TrajectoryStateOnSurface predTsos = thePropagator->propagate(tsosBefore, (**idet).surface());
00969 if (predTsos.isValid()) {
00970 GlobalVector propagationDistance = predTsos.globalPosition() - tsosBefore.globalPosition();
00971 if (propagationDistance.mag() > maxPropagationDistance) continue;
00972 dmmap.insert(std::make_pair(*idet, predTsos));
00973 }
00974 }
00975
00976 if(debug_) std::cout << "TRAJECTORY INTERSECTS " << dmmap.size() << " DETECTORS." << std::endl;
00977
00978
00979 std::map<const GeomDet*, TrajectoryMeasurement> dtmmap;
00980 for (edm::OwnVector<TrackingRecHit>::const_iterator ih = theHits.begin(); ih != theHits.end(); ++ih) {
00981 const GeomDet* det = trackerGeom->idToDet(ih->geographicalId());
00982
00983
00984
00985 std::map<const GeomDet*, TrajectoryStateOnSurface>::iterator idm = dmmap.find(det);
00986 if (idm == dmmap.end()) continue;
00987 TrajectoryStateOnSurface predTsos = idm->second;
00988 TransientTrackingRecHit::RecHitPointer rhit = ttrhBuilder->build(&(*ih));
00989
00990 const SiStripMatchedRecHit2D *origHit = dynamic_cast<const SiStripMatchedRecHit2D *>(&(*ih));
00991 if (origHit !=0){
00992 const GluedGeomDet *gdet = dynamic_cast<const GluedGeomDet*>(det);
00993 const SiStripMatchedRecHit2D *corrHit = theHitMatcher->match(origHit,gdet,predTsos.localDirection());
00994 if (corrHit!=0){
00995 rhit = ttrhBuilder->build(&(*corrHit));
00996 delete corrHit;
00997 }
00998 }
00999
01000 if (debug_) {
01001 std::cout << "predTSOS (x/y/z): " << predTsos.globalPosition().x()<< " / " << predTsos.globalPosition().y()<< " / " << predTsos.globalPosition().z() << std::endl;
01002 std::cout << "local position x: " << predTsos.localPosition().x() << "+-" << sqrt(predTsos.localError().positionError().xx()) << std::endl;
01003 std::cout << "local position y: " << predTsos.localPosition().y() << "+-" << sqrt(predTsos.localError().positionError().yy()) << std::endl;
01004 std::cout << "rhit local position x: " << rhit->localPosition().x() << "+-" << sqrt(rhit->localPositionError().xx()) << std::endl;
01005 std::cout << "rhit local position y: " << rhit->localPosition().y() << "+-" << sqrt(rhit->localPositionError().yy()) << std::endl;
01006 }
01007
01008 MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rhit);
01009 if (debug_) std::cout<< "hit " << ih-theHits.begin()
01010 << ": est = " << est.first << " " << est.second <<std::endl;
01011
01012
01013
01014 if (est.first) {
01015 TrajectoryMeasurement tm;
01016 TrajectoryStateOnSurface currTsos = theUpdator->update(predTsos, *rhit);
01017 std::map<const GeomDet*, TrajectoryMeasurement>::iterator idtm = dtmmap.find(det);
01018 if (idtm == dtmmap.end()) {
01019 tm = TrajectoryMeasurement (predTsos, currTsos, &(*rhit),est.second,
01020 theMeasurementTracker->geometricSearchTracker()->detLayer(ih->geographicalId()));
01021 dtmmap.insert(std::make_pair(det, tm));
01022 } else if (idtm->second.estimate() > est.second) {
01023 dtmmap.erase(idtm);
01024 tm = TrajectoryMeasurement(predTsos, currTsos, &(*rhit),est.second,
01025 theMeasurementTracker->geometricSearchTracker()->detLayer(ih->geographicalId()));
01026 dtmmap.insert(std::make_pair(det, tm));
01027 }
01028 if ((firstTM)){
01029 theBestTM = tm;
01030 if (debug_) std::cout <<"Initialize best to " << theBestTM.estimate() << std::endl;
01031 firstTM = false;
01032 }
01033 else if (!firstTM) {
01034 if (debug_) std::cout << "Current best is " << theBestTM.estimate() << " while this hit is " << est.second;
01035 if (est.second < theBestTM.estimate()) {
01036 if (debug_) std::cout << " so replace it " ;
01037 theBestTM = tm;
01038 }
01039 if (debug_) std::cout << std::endl;
01040 }
01041 }
01042 }
01043 if (debug_) std::cout<< "Hits(Dets) to add: " << dtmmap.size() <<std::endl;
01044 if (!dtmmap.empty()) {
01045
01046 std::vector<std::pair<TransientTrackingRecHit::ConstRecHitPointer, TrajectoryMeasurement*> > OverlapHits;
01047 for (std::map<const GeomDet*, TrajectoryMeasurement>::iterator idtm = dtmmap.begin();
01048 idtm != dtmmap.end(); ++idtm) {
01049 OverlapHits.push_back(std::make_pair(idtm->second.recHit(),&idtm->second));
01050
01051 if (debug_) std::cout<<" Measurement on layer "
01052 << theMeasurementTracker->geometricSearchTracker()->detLayer(idtm->second.recHit()->geographicalId())
01053 << " with estimate " << idtm->second.estimate()<<std::endl ;
01054 }
01055 if (debug_)
01056 std::cout<<" Best Measurement is on layer "
01057 << theMeasurementTracker->geometricSearchTracker()->detLayer(theBestTM.recHit()->geographicalId())
01058 << " with estimate " << theBestTM.estimate()<<std::endl ;
01059
01060
01061 if (dtmmap.size()==1){
01062 for (std::map<const GeomDet*, TrajectoryMeasurement>::iterator idtm = dtmmap.begin();
01063 idtm != dtmmap.end(); ++idtm) {
01064 TrajectoryMeasurement itm = idtm->second;
01065 if (debug_) std::cout<<" Measurement on layer "
01066 << theMeasurementTracker->geometricSearchTracker()->detLayer(itm.recHit()->geographicalId())
01067 << " with estimate " << itm.estimate()<<std::endl ;
01068 theBestHits.push_back(itm);
01069 }
01070 }
01071 else if (dtmmap.size()>=2) {
01072
01073 if (debug_) std::cout<<"Unsorted OverlapHits has size " <<OverlapHits.size() << std::endl;
01074
01075 for (std::vector<std::pair<TransientTrackingRecHit::ConstRecHitPointer,TrajectoryMeasurement*> >::iterator irh =OverlapHits.begin();
01076 irh!=OverlapHits.end();++irh){
01077 if (debug_) std::cout << "Hit " << irh-OverlapHits.begin()
01078 << " on det " << irh->first->det()
01079 << " detLayer "
01080 << theMeasurementTracker->geometricSearchTracker()->detLayer(irh->first->geographicalId())
01081 << ", r/phi/z = "
01082 << irh->first->globalPosition().perp() << " "
01083 << irh->first->globalPosition().phi() << " "
01084 << irh->first->globalPosition().z()
01085 << std::endl;
01086 }
01087
01088 std::sort( OverlapHits.begin(),OverlapHits.end(),SortHitTrajectoryPairsByGlobalPosition());
01089 if (debug_) std::cout<<"Sorted OverlapHits has size " <<OverlapHits.size() << std::endl;
01090
01091 float workingBestChi2 = 1000000.0;
01092 std::vector<TrajectoryMeasurement> workingBestHits;
01093
01094 std::vector<std::pair<TransientTrackingRecHit::ConstRecHitPointer,TrajectoryMeasurement*> >::iterator irh1;
01095 std::vector<std::pair<TransientTrackingRecHit::ConstRecHitPointer,TrajectoryMeasurement*> >::iterator irh2;
01096 for (irh1 =OverlapHits.begin(); irh1!=--OverlapHits.end(); ++irh1){
01097 theBestHits.clear();
01098 float running_chi2=0;
01099 if (debug_) std::cout << "Hit " << irh1-OverlapHits.begin()
01100 << " on det " << irh1->first->det()
01101 << " detLayer "
01102 << theMeasurementTracker->geometricSearchTracker()->detLayer(irh1->first->geographicalId())
01103 << ", r/phi/z = "
01104
01105 << irh1->first->globalPosition().perp() << " "
01106 << irh1->first->globalPosition().phi() << " "
01107 << irh1->first->globalPosition().z()
01108 << std::endl;
01109
01110 TrajectoryStateOnSurface currTsos = irh1->second->updatedState();
01111 TransientTrackingRecHit::ConstRecHitPointer rhit = irh1->first;
01112 theBestHits.push_back(*(irh1->second));
01113 if (debug_) std::cout<<"Added first hit with chi2 = " << irh1->second->estimate() << std::endl;
01114 running_chi2 += irh1->second->estimate();
01115 for (irh2 = irh1; irh2!=OverlapHits.end(); ++irh2){
01116 if (irh2 == irh1) continue;
01117 TransientTrackingRecHit::ConstRecHitPointer rh = irh2->first;
01118 const GeomDet* det = irh2->first->det();
01119
01120 TrajectoryStateOnSurface predTsos = thePropagator->propagate(currTsos, det->surface());
01121
01122 if (predTsos.isValid()){
01123 MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rh);
01124 if (debug_) std::cout<<"Added overlap hit with est = " << est.first << " " << est.second << std::endl;
01125 if (est.first){
01126 TrajectoryMeasurement tm(predTsos, currTsos, &(*rh),est.second,
01127 theMeasurementTracker->geometricSearchTracker()->detLayer(rh->geographicalId()));
01128 theBestHits.push_back(tm);
01129 running_chi2 += est.second ;
01130 }
01131 else {
01132 }
01133 }
01134
01135 }
01136 if (theBestHits.size()==dtmmap.size()){
01137 if (debug_) std::cout<<"Added all "<<theBestHits.size()<<" hits out of " << dtmmap.size() << std::endl;
01138 break;
01139 }
01140
01141 if (theBestHits.size() < dtmmap.size()){
01142 if (debug_) std::cout<<"Added only "<<theBestHits.size()<<" hits out of " << dtmmap.size() << std::endl;
01143
01144 if (theBestHits.size() > workingBestHits.size()){
01145 if (debug_) std::cout<<"Current combo has more hits so replace best" << std::endl;
01146 workingBestHits = theBestHits;
01147 }
01148
01149 else if (theBestHits.size() == workingBestHits.size()){
01150 if (running_chi2< workingBestChi2){
01151 if (debug_) std::cout<<"Current combo has same # of hits but lower chi2 so replace best" << std::endl;
01152 workingBestHits = theBestHits;
01153 workingBestChi2 = running_chi2;
01154 }
01155 }
01156 }
01157 }
01158 if (theBestHits.size()<2){
01159 if (debug_) std::cout<<"Only one good hit in overlap"<<std::endl;
01160 if (debug_) std::cout<<" Added hit on layer on det "
01161 << theBestTM.recHit()->det()
01162 << " detLayer "
01163 << theMeasurementTracker->geometricSearchTracker()->detLayer(theBestTM.recHit()->geographicalId())
01164 << ", r/phi/z = "
01165 << theBestTM.recHit()->globalPosition().perp() << " "
01166 << theBestTM.recHit()->globalPosition().phi() << " "
01167 << theBestTM.recHit()->globalPosition().z()
01168 << " with estimate " << theBestTM.estimate()<<std::endl ;
01169 theBestHits.clear();
01170 theBestHits.push_back(theBestTM);
01171 }
01172
01173 }
01174 else {
01175 if (debug_)std::cout << "ERROR: Unexpected size from DTMMAP = " << dtmmap.size() << std::endl;
01176 theBestHits.push_back(theBestTM);
01177 }
01178 }
01179
01180 return theBestHits;
01181 }
01182
01183
01184
01185 bool RoadSearchTrackCandidateMakerAlgorithm::chooseStartingLayers( std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >& recHitsByLayer,
01186
01187 std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilyr0,
01188 const std::multimap<int, const DetLayer*>& layer_map,
01189 std::set<const DetLayer*>& good_layers,
01190 std::vector<const DetLayer*>& middle_layers ,
01191 RoadSearchCloud::RecHitVector& recHits_middle)
01192 {
01193 const unsigned int max_middle_layers = 2;
01194
01195
01196
01197
01198
01199
01200 std::multimap<int, const DetLayer*>::const_iterator ilm = layer_map.begin();
01201 int ngoodlayers = 0;
01202 while (ilm != layer_map.end()) {
01203 if (ngoodlayers >= nFoundMin_ && ilm->first > 1) break;
01204
01205 good_layers.insert(ilm->second);
01206 ++ngoodlayers;
01207 ++ilm;
01208 }
01209
01210
01211 for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilayer = ilyr0+1;
01212 ilayer != recHitsByLayer.end(); ++ilayer) {
01213
01214 if (good_layers.find(ilayer->first) == good_layers.end()) continue;
01215
01216 if (!CosmicReco_ && !lstereo[ilayer-recHitsByLayer.begin()]) continue;
01217 middle_layers.push_back(ilayer->first);
01218 if (middle_layers.size() >= max_middle_layers) break;
01219 }
01220
01221 for (std::vector<const DetLayer*>::iterator ml = middle_layers.begin();
01222 ml!=middle_layers.end();++ml){
01223 unsigned int middle_layers_found = 0;
01224 for (std::vector<std::pair<const DetLayer*, RoadSearchCloud::RecHitVector > >::iterator ilyr = recHitsByLayer.begin();
01225 ilyr != recHitsByLayer.end(); ++ilyr) {
01226 if (ilyr->first == *ml){
01227 for (RoadSearchCloud::RecHitVector::const_iterator ih = ilyr->second.begin();
01228 ih != ilyr->second.end(); ++ih) {
01229 recHits_middle.push_back(*ih);
01230 }
01231 ++middle_layers_found;
01232 }
01233 if (middle_layers_found == middle_layers.size()) continue;
01234 }
01235
01236 }
01237 return (recHits_middle.size()>0);
01238 }
01239
01240 FreeTrajectoryState RoadSearchTrackCandidateMakerAlgorithm::initialTrajectory(const edm::EventSetup& es,
01241 const TrackingRecHit* theInnerHit,
01242 const TrackingRecHit* theOuterHit)
01243 {
01244 FreeTrajectoryState fts;
01245
01246 GlobalPoint inner = trackerGeom->idToDet(theInnerHit->geographicalId())->surface().toGlobal(theInnerHit->localPosition());
01247 GlobalPoint outer = trackerGeom->idToDet(theOuterHit->geographicalId())->surface().toGlobal(theOuterHit->localPosition());
01248
01249 LogDebug("RoadSearch") << "inner hit: r/phi/z = "<< inner.perp() << " " << inner.phi() << " " << inner.z() ;
01250 LogDebug("RoadSearch") << "outer hit: r/phi/z = "<< outer.perp() << " " << outer.phi() << " " << outer.z() ;
01251
01252
01253 const double dRmin = 0.1;
01254 if (outer.perp() - inner.perp() < dRmin) return fts;
01255
01256 const double dr2 = initialVertexErrorXY_*initialVertexErrorXY_;
01257
01258
01259 const double dz2 = 5.3*5.3;
01260
01261
01262 FastLine linearFit(outer, inner);
01263 double z_0 = -linearFit.c()/linearFit.n2();
01264 if ( std::abs(z_0) > 275 && !CosmicReco_ ) return fts;
01265
01266 GlobalError vertexErr(dr2,
01267 0, dr2,
01268 0, 0, dz2);
01269
01270
01271
01272
01273
01274 double x0=0.0,y0=0.0,z0=0.0;
01275 double phi0 = -999.0;
01276 if (NoFieldCosmic_){
01277 phi0=atan2(outer.y()-inner.y(),outer.x()-inner.x());
01278 thePropagator = theAloPropagator;
01279 if (inner.y()<outer.y()){
01280 if (debug_) std::cout<<"Flipping direction!!!" << std::endl;
01281 thePropagator = theRevPropagator;
01282 phi0=phi0-M_PI;
01283 }
01284 double alpha=atan2(inner.y(),inner.x());
01285 double d1=sqrt(inner.x()*inner.x()+inner.y()*inner.y());
01286 double d0=-d1*sin(alpha-phi0); x0=d0*sin(phi0); y0=-d0*cos(phi0);
01287 double l1=0.0,l2=0.0;
01288 if (fabs(cos(phi0))>0.1){
01289 l1=(inner.x()-x0)/cos(phi0);l2=(outer.x()-x0)/cos(phi0);
01290 }else{
01291 l1=(inner.y()-y0)/sin(phi0);l2=(outer.y()-y0)/sin(phi0);
01292 }
01293 z0=(l2*inner.z()-l1*outer.z())/(l2-l1);
01294 }
01295
01296 FastHelix helix(outer, inner, GlobalPoint(x0,y0,z0), es);
01297 if (!NoFieldCosmic_ && !helix.isValid()) return fts;
01298
01299 AlgebraicSymMatrix55 C = AlgebraicMatrixID();
01300 float zErr = vertexErr.czz();
01301 float transverseErr = vertexErr.cxx();
01302 C(3, 3) = transverseErr;
01303 C(4, 4) = zErr;
01304 CurvilinearTrajectoryError initialError(C);
01305
01306 if (NoFieldCosmic_) {
01307 TrackCharge q = 1;
01308 FastLine flfit(outer, inner);
01309 double dzdr = -flfit.n1()/flfit.n2();
01310 if (inner.y()<outer.y()) dzdr*=-1;
01311
01312 GlobalPoint XYZ0(x0,y0,z0);
01313 if (debug_) std::cout<<"Initial Point (x0/y0/z0): " << x0 <<'\t'<< y0 <<'\t'<< z0 << std::endl;
01314 GlobalVector PXYZ(cosmicSeedPt_*cos(phi0),
01315 cosmicSeedPt_*sin(phi0),
01316 cosmicSeedPt_*dzdr);
01317 GlobalTrajectoryParameters thePars(XYZ0,PXYZ,q,magField);
01318 AlgebraicSymMatrix66 CErr = AlgebraicMatrixID();
01319 CErr *= 5.0;
01320
01321
01322 fts = FreeTrajectoryState(thePars,
01323 CartesianTrajectoryError(CErr));
01324 if (debug_) std::cout<<"\nInitial CError (dx/dy/dz): " << CErr(1,1) <<'\t'<< CErr(2,2) <<'\t'<< CErr(3,3) << std::endl;
01325 if (debug_) std::cout<<"\n\ninner dy = " << theInnerHit->localPositionError().yy() <<"\t\touter dy = " << theOuterHit->localPositionError().yy() << std::endl;
01326 }
01327 else {
01328 fts = FreeTrajectoryState( helix.stateAtVertex().parameters(), initialError);
01329 }
01330
01331
01332
01333 return fts;
01334 }
01335
01336 FreeTrajectoryState RoadSearchTrackCandidateMakerAlgorithm::initialTrajectoryFromTriplet(const edm::EventSetup& es,
01337 const TrackingRecHit* theInnerHit,
01338 const TrackingRecHit* theMiddleHit,
01339 const TrackingRecHit* theOuterHit)
01340 {
01341 FreeTrajectoryState fts;
01342
01343 GlobalPoint inner = trackerGeom->idToDet(theInnerHit->geographicalId())->surface().toGlobal(theInnerHit->localPosition());
01344 GlobalPoint middle= trackerGeom->idToDet(theMiddleHit->geographicalId())->surface().toGlobal(theMiddleHit->localPosition());
01345 GlobalPoint outer = trackerGeom->idToDet(theOuterHit->geographicalId())->surface().toGlobal(theOuterHit->localPosition());
01346
01347 LogDebug("RoadSearch") << "inner hit: r/phi/z = "<< inner.perp() << " " << inner.phi() << " " << inner.z() ;
01348 LogDebug("RoadSearch") << "middlehit: r/phi/z = "<< inner.perp() << " " << inner.phi() << " " << inner.z() ;
01349 LogDebug("RoadSearch") << "outer hit: r/phi/z = "<< outer.perp() << " " << outer.phi() << " " << outer.z() ;
01350
01351
01352 const double dRmin = 0.1;
01353 if (outer.perp() - inner.perp() < dRmin) return fts;
01354 const double dr2 = initialVertexErrorXY_*initialVertexErrorXY_;
01355 const double dz2 = 5.3*5.3;
01356
01357
01358 FastLine linearFit(outer, inner);
01359 double z_0 = -linearFit.c()/linearFit.n2();
01360 if ( std::abs(z_0) > 275 && !CosmicReco_ ) return fts;
01361
01362
01363 FastHelix helix(outer, middle, inner, es);
01364
01365 if (!helix.isValid() ||
01366 std::isnan(helix.stateAtVertex().parameters().momentum().perp())) return fts;
01367
01368 if (helix.stateAtVertex().parameters().momentum().perp() < 0.5 ||
01369 std::abs(helix.stateAtVertex().parameters().position().z()) > 550.0)
01370 return fts;
01371
01372 AlgebraicSymMatrix55 C = AlgebraicMatrixID();
01373 float zErr = dz2;
01374 float transverseErr = dr2;
01375 C(3, 3) = transverseErr;
01376 C(4, 4) = zErr;
01377 CurvilinearTrajectoryError initialError(C);
01378
01379
01380 thePropagator = theAloPropagator;
01381 GlobalVector gv=helix.stateAtVertex().parameters().momentum();
01382 float charge=helix.stateAtVertex().parameters().charge();
01383
01384 if (CosmicReco_ && gv.y()>0){
01385 if (debug_) std::cout<<"Flipping direction!!!" << std::endl;
01386 thePropagator = theRevPropagator;
01387 gv=-1.*gv;
01388 charge=-1.*charge;
01389 }
01390
01391 GlobalTrajectoryParameters Gtp(inner,gv,int(charge),&(*magField));
01392 fts = FreeTrajectoryState(Gtp, initialError);
01393
01394
01395
01396 return fts;
01397 }
01398
01399
01400
01401 Trajectory RoadSearchTrackCandidateMakerAlgorithm::createSeedTrajectory(FreeTrajectoryState& fts,
01402 const TrackingRecHit* theInnerHit,
01403 const DetLayer* theInnerHitLayer)
01404
01405 {
01406 Trajectory theSeedTrajectory;
01407
01408
01409 const GeomDet* innerDet = trackerGeom->idToDet((theInnerHit)->geographicalId());
01410 const TrajectoryStateOnSurface innerState =
01411 thePropagator->propagate(fts,innerDet->surface());
01412 if ( !innerState.isValid() ||
01413 std::isnan(innerState.globalMomentum().perp())) {
01414 if (debug_) std::cout<<"*******DISASTER ********* seed doesn't make it to first hit!!!!!" << std::endl;
01415 return theSeedTrajectory;
01416 }
01417 TransientTrackingRecHit::RecHitPointer intrhit = ttrhBuilder->build(theInnerHit);
01418
01419 const SiStripMatchedRecHit2D *origHit = dynamic_cast<const SiStripMatchedRecHit2D *>(theInnerHit);
01420 if (origHit !=0){
01421 const GluedGeomDet *gdet = dynamic_cast<const GluedGeomDet*>(innerDet);
01422 const SiStripMatchedRecHit2D *corrHit = theHitMatcher->match(origHit,gdet,innerState.localDirection());
01423 if (corrHit!=0){
01424 intrhit = ttrhBuilder->build(&(*corrHit));
01425 delete corrHit;
01426 }
01427 }
01428
01429 MeasurementEstimator::HitReturnType est = theEstimator->estimate(innerState, *intrhit);
01430 if (!est.first) return theSeedTrajectory;
01431 TrajectoryStateOnSurface innerUpdated= theUpdator->update( innerState,*intrhit);
01432 if (debug_) std::cout<<"InnerUpdated: " << innerUpdated << std::endl;
01433 if (!innerUpdated.isValid() ||
01434 std::isnan(innerUpdated.globalMomentum().perp())) {
01435 if (debug_) std::cout<<"Trajectory updated with first hit is invalid!!!" << std::endl;
01436 return theSeedTrajectory;
01437 }
01438 TrajectoryMeasurement tm = TrajectoryMeasurement(innerState, innerUpdated, &(*intrhit),est.second,theInnerHitLayer);
01439
01440 PTrajectoryStateOnDet pFirstStateTwo = trajectoryStateTransform::persistentState(innerUpdated,
01441 intrhit->geographicalId().rawId());
01442 edm::OwnVector<TrackingRecHit> newHitsTwo;
01443 newHitsTwo.push_back(intrhit->hit()->clone());
01444
01445 TrajectorySeed tmpseedTwo = TrajectorySeed(pFirstStateTwo,
01446 newHitsTwo,
01447 alongMomentum);
01448 if (thePropagator->propagationDirection()==oppositeToMomentum) {
01449 tmpseedTwo = TrajectorySeed(pFirstStateTwo,
01450 newHitsTwo,
01451 oppositeToMomentum);
01452 }
01453
01454
01455
01456 theSeedTrajectory = Trajectory(tmpseedTwo, tmpseedTwo.direction());
01457
01458 theSeedTrajectory.push(tm,est.second);
01459
01460 return theSeedTrajectory;
01461 }
01462
01463
01464
01465 std::vector<Trajectory> RoadSearchTrackCandidateMakerAlgorithm::extrapolateTrajectory(const Trajectory& theTrajectory,
01466 RoadSearchCloud::RecHitVector& theLayerHits,
01467 const DetLayer* innerHitLayer,
01468 const TrackingRecHit* outerHit,
01469 const DetLayer* outerHitLayer)
01470 {
01471 std::vector<Trajectory> newTrajectories;
01472
01473 for(RoadSearchCloud::RecHitVector::const_iterator ihit = theLayerHits.begin();
01474 ihit != theLayerHits.end(); ihit++) {
01475 Trajectory traj = theTrajectory;
01476 const DetLayer* thisLayer =
01477 theMeasurementTracker->geometricSearchTracker()->detLayer((*ihit)->geographicalId());
01478 if (thisLayer == innerHitLayer){
01479
01480
01481
01482 continue;
01483 }
01484 if (thisLayer == outerHitLayer){
01485 LocalPoint p1 = (*ihit)->localPosition();
01486 LocalPoint p2 = outerHit->localPosition();
01487 if (p1.x()!=p2.x() || p1.y()!=p2.y()) continue;
01488 }
01489
01490
01491 TransientTrackingRecHit::RecHitPointer rhit = ttrhBuilder->build(*ihit);
01492
01493 if (debug_){
01494 if (rhit->isValid()) {
01495 LogDebug("RoadSearch") << "RecHit " << ihit-theLayerHits.begin()
01496 << ", det " << rhit->det() << ", r/phi/z = "
01497 << rhit->globalPosition().perp() << " "
01498 << rhit->globalPosition().phi() << " "
01499 << rhit->globalPosition().z();
01500 } else {
01501 LogDebug("RoadSearch") << "RecHit " << ihit-theLayerHits.begin()
01502 << " (invalid)";
01503 }
01504 }
01505
01506 const GeomDet* det = trackerGeom->idToDet(rhit->geographicalId());
01507
01508 TrajectoryStateOnSurface predTsos;
01509 TrajectoryStateOnSurface currTsos;
01510
01511 if (traj.measurements().empty()) {
01512
01513 if (debug_) std::cout<<"BIG ERROR!!! How did we make it to here with no trajectory measurements?!?!?"<<std::endl;
01514 } else {
01515 currTsos = traj.measurements().back().updatedState();
01516 predTsos = thePropagator->propagate(currTsos, det->surface());
01517 }
01518 if (!predTsos.isValid()){
01519 continue;
01520 }
01521 GlobalVector propagationDistance = predTsos.globalPosition() - currTsos.globalPosition();
01522 if (propagationDistance.mag() > maxPropagationDistance) continue;
01523 if (debug_) {
01524 std::cout << "currTsos (x/y/z): "
01525 << currTsos.globalPosition().x() << " / "
01526 << currTsos.globalPosition().y() << " / "
01527 << currTsos.globalPosition().z() << std::endl;
01528 std::cout << "predTsos (x/y/z): "
01529 << predTsos.globalPosition().x() << " / "
01530 << predTsos.globalPosition().y() << " / "
01531 << predTsos.globalPosition().z() << std::endl;
01532 std::cout << "Propagation distance1 is " << propagationDistance.mag() << std::endl;
01533 }
01534 TrajectoryMeasurement tm;
01535 if (debug_){
01536 std::cout<< "trajectory at r/z=" << det->surface().position().perp()
01537 << " " << det->surface().position().z()
01538 << ", hit " << ihit-theLayerHits.begin()
01539 << " local prediction " << predTsos.localPosition().x()
01540 << " +- " << sqrt(predTsos.localError().positionError().xx())
01541 << ", hit at " << rhit->localPosition().x() << " +- " << sqrt(rhit->localPositionError().xx())
01542 << std::endl;
01543 }
01544
01545
01546
01547
01548 const SiStripMatchedRecHit2D *origHit = dynamic_cast<const SiStripMatchedRecHit2D *>(*ihit);
01549 if (origHit !=0){
01550 const GluedGeomDet *gdet = dynamic_cast<const GluedGeomDet*>(rhit->det());
01551 const SiStripMatchedRecHit2D *corrHit = theHitMatcher->match(origHit,gdet,predTsos.localDirection());
01552 if (corrHit!=0){
01553 rhit = ttrhBuilder->build(&(*corrHit));
01554 delete corrHit;
01555 }
01556 }
01557
01558 MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rhit);
01559 if (debug_) std::cout << "estimation: " << est.first << " " << est.second << std::endl;
01560 if (!est.first) continue;
01561 currTsos = theUpdator->update(predTsos, *rhit);
01562 tm = TrajectoryMeasurement(predTsos, currTsos, &(*rhit),est.second,thisLayer);
01563 traj.push(tm,est.second);
01564 newTrajectories.push_back(traj);
01565
01566
01567 }
01568
01569 return newTrajectories;
01570 }
01571
01572
01573 TrackCandidateCollection RoadSearchTrackCandidateMakerAlgorithm::PrepareTrackCandidates(std::vector<Trajectory>& theTrajectories)
01574 {
01575
01576 TrackCandidateCollection theCollection;
01577
01578 theTrajectoryCleaner->clean(theTrajectories);
01579
01580
01581
01582 if(CosmicTrackMerging_) {
01583
01584
01585 std::vector<Trajectory> traj;
01586
01587
01588 std::vector<bool> trajUsed;
01589
01590 for (std::vector<Trajectory>::iterator it = theTrajectories.begin(); it != theTrajectories.end(); ++it) {
01591 if (it->isValid()) {
01592 traj.push_back(*it);
01593 trajUsed.push_back(false);
01594 }
01595 }
01596
01597 if(debugCosmics_) {
01598 std::cout << "==========ENTERING COSMIC MODE===========" << std::endl;
01599
01600 for (std::vector<Trajectory>::iterator it = traj.begin(); it != traj.end(); it++) {
01601 std::cout << "Trajectory " << it-traj.begin() << " has "<<it->recHits().size()<<" hits and is valid: " << it->isValid() << std::endl;
01602 TransientTrackingRecHit::ConstRecHitContainer itHits = it->recHits();
01603
01604
01605 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator rhit=itHits.begin(); rhit!=itHits.end(); ++rhit)
01606 std::cout << "-->good hit position: " << (*rhit)->globalPosition().x() << ", "
01607 << (*rhit)->globalPosition().y() << ", "<< (*rhit)->globalPosition().z() << std::endl;
01608
01609 }
01610 }
01611
01612
01613 for ( unsigned int i = 0; i < traj.size(); ++i) {
01614 if (trajUsed[i]) continue;
01615 for ( unsigned int j = i+1; j != traj.size(); ++j) {
01616 if (trajUsed[j]) continue;
01617
01618 if (debugCosmics_) std::cout<< "Trajectory " <<i<< " has "<<traj[i].recHits().size()<<" hits with chi2=" << traj[i].chiSquared() << " and is valid"<<std::endl;
01619 if (debugCosmics_) std::cout<< "Trajectory " <<j<< " has "<<traj[j].recHits().size()<<" hits with chi2=" << traj[j].chiSquared() << " and is valid"<<std::endl;
01620
01621 TrajectoryMeasurement firstTraj1 = traj[i].firstMeasurement();
01622 TrajectoryMeasurement firstTraj2 = traj[j].firstMeasurement();
01623 TrajectoryStateOnSurface firstTraj1TSOS = firstTraj1.updatedState();
01624 TrajectoryStateOnSurface firstTraj2TSOS = firstTraj2.updatedState();
01625
01626
01627 if(debugCosmics_)
01628 std::cout << "phi1: " << firstTraj1TSOS.globalMomentum().phi()
01629 << " phi2: " << firstTraj2TSOS.globalMomentum().phi()
01630 << " --> delta_phi: " << firstTraj1TSOS.globalMomentum().phi()-firstTraj2TSOS.globalMomentum().phi() << std::endl;
01631
01632
01633
01634 if( fabs(firstTraj1TSOS.globalMomentum().phi()-firstTraj2TSOS.globalMomentum().phi())<0.3 ) {
01635 if(debugCosmics_) std::cout << "-->match successful" << std::endl;
01636 } else {
01637 if(debugCosmics_) std::cout << "-->match not successful" << std::endl;
01638 }
01639 if( fabs(firstTraj1TSOS.globalMomentum().phi()-firstTraj2TSOS.globalMomentum().phi())>0.3 ) continue;
01640
01641
01642
01643
01644 TrajectoryMeasurement lastTraj1 = traj[i].lastMeasurement();
01645 TrajectoryMeasurement lastTraj2 = traj[j].lastMeasurement();
01646 TrajectoryStateOnSurface lastTraj1TSOS = lastTraj1.updatedState();
01647 TrajectoryStateOnSurface lastTraj2TSOS = lastTraj2.updatedState();
01648
01649 if(debugCosmics_){
01650 std::cout<<"Traj1 first (x/y/z): "
01651 << firstTraj1TSOS.globalPosition().x() <<" / "
01652 << firstTraj1TSOS.globalPosition().y() <<" / "
01653 << firstTraj1TSOS.globalPosition().z()
01654 << " phi: " << firstTraj1TSOS.globalMomentum().phi() << std::endl;
01655 std::cout<<"Traj1 last (x/y/z): "
01656 << lastTraj1TSOS.globalPosition().x() <<" / "
01657 << lastTraj1TSOS.globalPosition().y() <<" / "
01658 << lastTraj1TSOS.globalPosition().z()
01659 << " phi: " << lastTraj1TSOS.globalMomentum().phi() << std::endl;
01660
01661 std::cout<<"Traj2 first (x/y/z): "
01662 << firstTraj2TSOS.globalPosition().x() <<" / "
01663 << firstTraj2TSOS.globalPosition().y() <<" / "
01664 << firstTraj2TSOS.globalPosition().z()
01665 << " phi: " << firstTraj2TSOS.globalMomentum().phi() << std::endl;
01666 std::cout<<"Traj2 last (x/y/z): "
01667 << lastTraj2TSOS.globalPosition().x() <<" / "
01668 << lastTraj2TSOS.globalPosition().y() <<" / "
01669 << lastTraj2TSOS.globalPosition().z()
01670 << " phi: " << lastTraj2TSOS.globalMomentum().phi() << std::endl;
01671
01672 }
01673
01674 Trajectory *upperTrajectory, *lowerTrajectory;
01675
01676 TrajectoryStateOnSurface lowerTSOS1,upperTSOS1;
01677 if (lastTraj1TSOS.globalPosition().y()<firstTraj1TSOS.globalPosition().y()) {
01678 lowerTSOS1=lastTraj1TSOS; upperTSOS1=firstTraj1TSOS;
01679 }
01680 else {
01681 lowerTSOS1=firstTraj1TSOS; upperTSOS1=lastTraj1TSOS;
01682 }
01683 TrajectoryStateOnSurface lowerTSOS2;
01684 if (lastTraj2TSOS.globalPosition().y()<firstTraj2TSOS.globalPosition().y()) lowerTSOS2=lastTraj2TSOS;
01685 else lowerTSOS2=firstTraj2TSOS;
01686 if(lowerTSOS1.globalPosition().y() > lowerTSOS2.globalPosition().y()) {
01687 if(debugCosmics_)
01688 std::cout << "-->case A: "<< lowerTSOS1.globalPosition().y() << " > " << lowerTSOS2.globalPosition().y() << std::endl;
01689
01690 upperTrajectory = &(traj[i]);
01691 lowerTrajectory = &(traj[j]);
01692
01693 } else {
01694 if(debugCosmics_)
01695 std::cout << "-->case B: "<< lowerTSOS1.globalPosition().y() << " < " << lowerTSOS2.globalPosition().y() << std::endl;
01696
01697 upperTrajectory = &(traj[j]);
01698 lowerTrajectory = &(traj[i]);
01699 }
01700
01701 std::vector<Trajectory> freshStartUpperTrajectory = theSmoother->trajectories(*upperTrajectory);
01702 std::vector<Trajectory> freshStartLowerTrajectory = theSmoother->trajectories(*lowerTrajectory);
01703
01704 if (freshStartUpperTrajectory.empty() || freshStartLowerTrajectory .empty()){
01705 if (debugCosmics_) std::cout << " the smoother has failed."<<std::endl;
01706 continue;
01707 }
01708
01709 TrajectoryStateOnSurface NewFirstTsos = freshStartUpperTrajectory.begin()->lastMeasurement().updatedState();
01710 TrajectoryStateOnSurface forwardTsos = freshStartUpperTrajectory.begin()->firstMeasurement().forwardPredictedState();
01711 TrajectoryStateOnSurface backwardTsos = freshStartUpperTrajectory.begin()->lastMeasurement().backwardPredictedState();
01712
01713 Trajectory freshStartTrajectory = *freshStartUpperTrajectory.begin();
01714 if(debugCosmics_) std::cout << "seed hit updatedState: " << NewFirstTsos.globalMomentum().x() << ", " << NewFirstTsos.globalMomentum().y() << ", " << NewFirstTsos.globalMomentum().z() << std::endl;
01715 if(debugCosmics_) std::cout << "seed hit updatedState (pos x/y/z): " << NewFirstTsos.globalPosition().x() << ", " << NewFirstTsos.globalPosition().y() << ", " << NewFirstTsos.globalPosition().z() << std::endl;
01716 if(debugCosmics_) std::cout << "seed hit forwardPredictedState: " << forwardTsos.globalMomentum().x() << ", " << forwardTsos.globalMomentum().y() << ", " << forwardTsos.globalMomentum().z() << std::endl;
01717 if(debugCosmics_) std::cout << "seed hit forwardPredictedState (pos x/y/z): " << forwardTsos.globalPosition().x() << ", " << forwardTsos.globalPosition().y() << ", " << forwardTsos.globalPosition().z() << std::endl;
01718 if(debugCosmics_) std::cout << "seed hit backwardPredictedState: " << backwardTsos.globalMomentum().x() << ", " << backwardTsos.globalMomentum().y() << ", " << backwardTsos.globalMomentum().z() << std::endl;
01719 if(debugCosmics_) std::cout << "seed hit backwardPredictedState (pos x/y/z): " << backwardTsos.globalPosition().x() << ", " << backwardTsos.globalPosition().y() << ", " << backwardTsos.globalPosition().z() << std::endl;
01720
01721 if(debugCosmics_) std::cout<<"#hits for upper trajectory: " << freshStartTrajectory.measurements().size() << std::endl;
01722
01723
01724 TransientTrackingRecHit::ConstRecHitContainer ttHits = lowerTrajectory->recHits(splitMatchedHits_);
01725
01726 if(debugCosmics_) std::cout << "loop over hits in lower trajectory..." << std::endl;
01727
01728 bool addHitToFreshStartTrajectory = false;
01729 bool propagationFailed = false;
01730 int lostHits = 0;
01731 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator rhit=ttHits.begin(); rhit!=ttHits.end(); ++rhit){
01732
01733 if(debugCosmics_ && lostHits>0){
01734 std::cout << " Lost " << lostHits << " of " << ttHits.size() << " on lower trajectory " << std::endl;
01735 std::cout << " Lost " << ((float)lostHits/(float)ttHits.size()) << " of hits of on lower trajectory " << std::endl;
01736 }
01737 if ((float)lostHits/(float)ttHits.size() > 0.5) {
01738 propagationFailed = true;
01739 break;
01740 }
01741 if(debugCosmics_) std::cout << "-->hit position: " << (*rhit)->globalPosition().x() << ", " << (*rhit)->globalPosition().y() << ", "<< (*rhit)->globalPosition().z() << std::endl;
01742
01743 TrajectoryStateOnSurface predTsos;
01744 TrajectoryStateOnSurface currTsos;
01745 TrajectoryMeasurement tm;
01746
01747 TransientTrackingRecHit::ConstRecHitPointer rh = (*rhit);
01748
01749
01750
01751
01752
01753
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768
01769
01770
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811 if( rh->isValid() ) {
01812 if(debugCosmics_) std::cout << "-->hit is valid"<<std::endl;
01813 const GeomDet* det = trackerGeom->idToDet(rh->geographicalId());
01814 if( addHitToFreshStartTrajectory==false ) {
01815
01816 currTsos = freshStartTrajectory.lastMeasurement().backwardPredictedState();
01817 } else {
01818
01819 currTsos = freshStartTrajectory.lastMeasurement().forwardPredictedState();
01820 }
01821
01822 if(debugCosmics_) std::cout << "current TSOS: " << currTsos.globalPosition().x() << ", " << currTsos.globalPosition().y() << ", " << currTsos.globalPosition().z() << ", " << std::endl;
01823
01824 predTsos = theAloPropagator->propagate(currTsos, det->surface());
01825
01826 if (!predTsos.isValid()) {
01827 if(debugCosmics_) std::cout<<"predTsos is not valid!" <<std::endl;
01828
01829 ++lostHits;
01830
01831 continue;
01832 }
01833 GlobalVector propagationDistance = predTsos.globalPosition() - currTsos.globalPosition();
01834 if (propagationDistance.mag() > maxPropagationDistance) continue;
01835
01836 if(debugCosmics_) std::cout << "predicted TSOS: " << predTsos.globalPosition().x() << ", " << predTsos.globalPosition().y() << ", " << predTsos.globalPosition().z() << ", " << std::endl;
01837 MeasurementEstimator::HitReturnType est = theEstimator->estimate(predTsos, *rh);
01838 if (!est.first) {
01839 if(debugCosmics_) std::cout<<"Failed to add one of the original hits on a low occupancy layer!!!!" << std::endl;
01840
01841 ++lostHits;
01842
01843 continue;
01844 }
01845
01846 currTsos = theUpdator->update(predTsos, *rh);
01847 if(debugCosmics_) std::cout << "current updated TSOS: " << currTsos.globalPosition().x() << ", " << currTsos.globalPosition().y() << ", " << currTsos.globalPosition().z() << ", " << std::endl;
01848 tm = TrajectoryMeasurement(predTsos, currTsos,&(*rh),est.second,theMeasurementTracker->geometricSearchTracker()->detLayer(rh->geographicalId()));
01849 freshStartTrajectory.push(tm,est.second);
01850 addHitToFreshStartTrajectory=true;
01851
01852 }
01853
01854 if(debugCosmics_) std::cout<<"#hits for new trajectory (his from upper trajectory added): " << freshStartTrajectory.measurements().size() << std::endl;
01855 }
01856
01857 if (propagationFailed) {
01858 if (debugCosmics_) std::cout<<"Propagation failed so go to next trajectory" << std::endl;
01859 continue;
01860 }
01861
01862
01863 if(debugCosmics_) std::cout << "put final trajectory together..." << std::endl;
01864 edm::OwnVector<TrackingRecHit> goodHits;
01865 TransientTrackingRecHit::ConstRecHitContainer tttempHits = freshStartTrajectory.recHits(splitMatchedHits_);
01866
01867 for (int k=tttempHits.size()-1; k>=0; k--) {
01868 if(debugCosmics_) std::cout << "-->good hit position: " << tttempHits[k]->globalPosition().x() << ", " << tttempHits[k]->globalPosition().y() << ", "<< tttempHits[k]->globalPosition().z() << std::endl;
01869 goodHits.push_back(tttempHits[k]->hit()->clone());
01870 }
01871 TrajectoryStateOnSurface firstState;
01872
01873
01874
01875 DetId FirstHitId = (*(freshStartTrajectory.recHits().end()-1))->geographicalId();
01876
01877 TrajectoryMeasurement maxYMeasurement = freshStartTrajectory.lastMeasurement();
01878 const GeomDet* det = trackerGeom->idToDet(FirstHitId);
01879 firstState = theAnalyticalPropagator->propagate(maxYMeasurement.updatedState(),det->surface());
01880 if (firstState.isValid() == false) continue;
01881 PTrajectoryStateOnDet const & state = trajectoryStateTransform::persistentState(firstState,FirstHitId.rawId());
01882
01883
01884 TrajectoryStateOnSurface firstTSOS = freshStartTrajectory.lastMeasurement().updatedState();
01885 if(debugCosmics_) std::cout << "generate new trajectory seed with hit (x/y/z): " << firstTSOS.globalPosition().x() << ", " << firstTSOS.globalPosition().y() << ", " << firstTSOS.globalPosition().z() << ", " << std::endl;
01886 TransientTrackingRecHit::ConstRecHitPointer rhit = freshStartTrajectory.lastMeasurement().recHit();
01887 PTrajectoryStateOnDet const & pFirstState = trajectoryStateTransform::persistentState(NewFirstTsos,rhit->geographicalId().rawId());
01888 edm::OwnVector<TrackingRecHit> newHits;
01889 newHits.push_back(rhit->hit()->clone());
01890 TrajectorySeed tmpseed = TrajectorySeed(pFirstState,newHits,alongMomentum);
01891
01892 theCollection.push_back(TrackCandidate(goodHits,freshStartTrajectory.seed(),state));
01893
01894
01895 trajUsed[i]=true;
01896 trajUsed[j]=true;
01897
01898 }
01899
01900 }
01901
01902
01903 for ( unsigned int i = 0; i < traj.size(); ++i) {
01904
01905 if (trajUsed[i]==true) continue;
01906
01907 if (debugCosmics_) std::cout<< "Trajectory (not merged) has "<<traj[i].recHits().size()<<" hits with chi2=" << traj[i].chiSquared() << " and is valid? "<< traj[i].isValid()<<std::endl;
01908 edm::OwnVector<TrackingRecHit> goodHits;
01909 TransientTrackingRecHit::ConstRecHitContainer ttHits = traj[i].recHits(splitMatchedHits_);
01910 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator rhit=ttHits.begin(); rhit!=ttHits.end(); ++rhit){
01911 goodHits.push_back((*rhit)->hit()->clone());
01912 }
01913 TrajectoryStateOnSurface firstState;
01914
01915
01916
01917 DetId FirstHitId = (*(traj[i].recHits().begin()))->geographicalId();
01918
01919
01920 TrajectoryMeasurement firstMeasurement = traj[i].measurements()[0];
01921 const GeomDet* det = trackerGeom->idToDet(FirstHitId);
01922 firstState = theAnalyticalPropagator->propagate(firstMeasurement.updatedState(), det->surface());
01923 if (firstState.isValid()) {
01924 PTrajectoryStateOnDet state = trajectoryStateTransform::persistentState(firstState,FirstHitId.rawId());
01925 theCollection.push_back(TrackCandidate(goodHits,traj[i].seed(),state));
01926 }
01927 }
01928 if (debugCosmics_) std::cout << "Original collection had " << theTrajectories.size()
01929 << " candidates. Merged collection has " << theCollection.size() << std::endl;
01930 }
01931
01932
01933 if(!CosmicTrackMerging_) {
01934 for (std::vector<Trajectory>::const_iterator it = theTrajectories.begin(); it != theTrajectories.end(); it++) {
01935 if (debug_) std::cout<< " Trajectory has "<<it->recHits().size()<<" hits with chi2=" << it->chiSquared() << " and is valid? "<<it->isValid()<<std::endl;
01936 if (it->isValid()){
01937
01938 edm::OwnVector<TrackingRecHit> goodHits;
01939 TransientTrackingRecHit::ConstRecHitContainer ttHits = it->recHits(splitMatchedHits_);
01940 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator rhit=ttHits.begin(); rhit!=ttHits.end(); ++rhit){
01941 goodHits.push_back((*rhit)->hit()->clone());
01942 }
01943 TrajectoryStateOnSurface firstState;
01944
01945
01946
01947 DetId FirstHitId = (*(it->recHits().begin()))->geographicalId();
01948
01949
01950 TrajectoryMeasurement firstMeasurement = it->measurements()[0];
01951 const GeomDet* det = trackerGeom->idToDet(FirstHitId);
01952 firstState = theAnalyticalPropagator->propagate(firstMeasurement.updatedState(), det->surface());
01953 if (firstState.isValid() == false) continue;
01954 PTrajectoryStateOnDet state = trajectoryStateTransform::persistentState(firstState,FirstHitId.rawId());
01955 theCollection.push_back(TrackCandidate(goodHits,it->seed(),state));
01956 }
01957 }
01958 }
01959
01960 return theCollection;
01961 }
1.7.3