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