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