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