00001 #include <memory>
00002
00003 #include "FWCore/Framework/interface/Event.h"
00004 #include "FWCore/Framework/interface/EventSetup.h"
00005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00006 #include "FWCore/Framework/interface/ESHandle.h"
00007
00008 #include "DataFormats/Common/interface/Handle.h"
00009 #include "DataFormats/Common/interface/OwnVector.h"
00010 #include "DataFormats/TrackerRecHit2D/interface/SiTrackerGSRecHit2DCollection.h"
00011 #include "DataFormats/TrackerRecHit2D/interface/SiTrackerGSMatchedRecHit2DCollection.h"
00012 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
00013 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00014 #include "DataFormats/VertexReco/interface/Vertex.h"
00015
00016 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00017 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00018
00019 #include "FastSimulation/ParticlePropagator/interface/MagneticFieldMapRecord.h"
00020
00021 #include "FastSimulation/Tracking/plugins/TrajectorySeedProducer.h"
00022 #include "FastSimulation/Tracking/interface/TrackerRecHit.h"
00023
00024 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00025 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
00026
00027 #include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
00028 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00029 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00030
00031 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
00032 #include "DataFormats/DetId/interface/DetId.h"
00033
00034 #include "FastSimulation/BaseParticlePropagator/interface/BaseParticlePropagator.h"
00035 #include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
00036
00037
00038 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048 TrajectorySeedProducer::TrajectorySeedProducer(const edm::ParameterSet& conf) :thePropagator(0)
00049 {
00050
00051
00052 theBeamSpot = conf.getParameter<edm::InputTag>("beamSpot");
00053
00054
00055 seedingAlgo = conf.getParameter<std::vector<std::string> >("seedingAlgo");
00056 for ( unsigned i=0; i<seedingAlgo.size(); ++i )
00057 produces<TrajectorySeedCollection>(seedingAlgo[i]);
00058
00059
00060 pTMin = conf.getParameter<std::vector<double> >("pTMin");
00061 if ( pTMin.size() != seedingAlgo.size() )
00062 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00063 << " WARNING : pTMin does not have the proper size "
00064 << std::endl;
00065
00066 for ( unsigned i=0; i<pTMin.size(); ++i )
00067 pTMin[i] *= pTMin[i];
00068
00069
00070 minRecHits = conf.getParameter<std::vector<unsigned int> >("minRecHits");
00071 if ( minRecHits.size() != seedingAlgo.size() )
00072 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00073 << " WARNING : minRecHits does not have the proper size "
00074 << std::endl;
00075
00076 absMinRecHits = 0;
00077 for ( unsigned ialgo=0; ialgo<minRecHits.size(); ++ialgo )
00078 if ( minRecHits[ialgo] > absMinRecHits ) absMinRecHits = minRecHits[ialgo];
00079
00080
00081 maxD0 = conf.getParameter<std::vector<double> >("maxD0");
00082 if ( maxD0.size() != seedingAlgo.size() )
00083 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00084 << " WARNING : maxD0 does not have the proper size "
00085 << std::endl;
00086
00087 maxZ0 = conf.getParameter<std::vector<double> >("maxZ0");
00088 if ( maxZ0.size() != seedingAlgo.size() )
00089 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00090 << " WARNING : maxZ0 does not have the proper size "
00091 << std::endl;
00092
00093
00094 hitProducer = conf.getParameter<edm::InputTag>("HitProducer");
00095
00096
00097 seedCleaning = conf.getParameter<bool>("seedCleaning");
00098
00099
00100 numberOfHits = conf.getParameter<std::vector<unsigned int> >("numberOfHits");
00101 if ( numberOfHits.size() != seedingAlgo.size() )
00102 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00103 << " WARNING : numberOfHits does not have the proper size "
00104 << std::endl;
00105
00106
00107 selectMuons = conf.getParameter<bool>("selectMuons");
00108
00109
00110 newSyntax = conf.getParameter<bool>("newSyntax");
00111 if (newSyntax) {
00112 layerList = conf.getParameter<std::vector<std::string> >("layerList");
00113
00114 } else {
00115
00116 firstHitSubDetectorNumber =
00117 conf.getParameter<std::vector<unsigned int> >("firstHitSubDetectorNumber");
00118 if ( firstHitSubDetectorNumber.size() != seedingAlgo.size() )
00119 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00120 << " WARNING : firstHitSubDetectorNumber does not have the proper size "
00121 << std::endl;
00122
00123 std::vector<unsigned int> firstSubDets =
00124 conf.getParameter<std::vector<unsigned int> >("firstHitSubDetectors");
00125 unsigned isub1 = 0;
00126 unsigned check1 = 0;
00127 firstHitSubDetectors.resize(seedingAlgo.size());
00128 for ( unsigned ialgo=0; ialgo<firstHitSubDetectorNumber.size(); ++ialgo ) {
00129 check1 += firstHitSubDetectorNumber[ialgo];
00130 for ( unsigned idet=0; idet<firstHitSubDetectorNumber[ialgo]; ++idet ) {
00131 firstHitSubDetectors[ialgo].push_back(firstSubDets[isub1++]);
00132 }
00133 }
00134 if ( firstSubDets.size() != check1 )
00135 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00136 << " WARNING : firstHitSubDetectors does not have the proper size (should be " << check1 << ")"
00137 << std::endl;
00138
00139
00140 secondHitSubDetectorNumber =
00141 conf.getParameter<std::vector<unsigned int> >("secondHitSubDetectorNumber");
00142 if ( secondHitSubDetectorNumber.size() != seedingAlgo.size() )
00143 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00144 << " WARNING : secondHitSubDetectorNumber does not have the proper size "
00145 << std::endl;
00146
00147 std::vector<unsigned int> secondSubDets =
00148 conf.getParameter<std::vector<unsigned int> >("secondHitSubDetectors");
00149 unsigned isub2 = 0;
00150 unsigned check2 = 0;
00151 secondHitSubDetectors.resize(seedingAlgo.size());
00152 for ( unsigned ialgo=0; ialgo<secondHitSubDetectorNumber.size(); ++ialgo ) {
00153 check2 += secondHitSubDetectorNumber[ialgo];
00154 for ( unsigned idet=0; idet<secondHitSubDetectorNumber[ialgo]; ++idet ) {
00155 secondHitSubDetectors[ialgo].push_back(secondSubDets[isub2++]);
00156 }
00157 }
00158 if ( secondSubDets.size() != check2 )
00159 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00160 << " WARNING : secondHitSubDetectors does not have the proper size (should be " << check2 << ")"
00161 << std::endl;
00162
00163 thirdHitSubDetectorNumber =
00164 conf.getParameter<std::vector<unsigned int> >("thirdHitSubDetectorNumber");
00165 if ( thirdHitSubDetectorNumber.size() != seedingAlgo.size() )
00166 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00167 << " WARNING : thirdHitSubDetectorNumber does not have the proper size "
00168 << std::endl;
00169
00170 std::vector<unsigned int> thirdSubDets =
00171 conf.getParameter<std::vector<unsigned int> >("thirdHitSubDetectors");
00172 unsigned isub3 = 0;
00173 unsigned check3 = 0;
00174 thirdHitSubDetectors.resize(seedingAlgo.size());
00175 for ( unsigned ialgo=0; ialgo<thirdHitSubDetectorNumber.size(); ++ialgo ) {
00176 check3 += thirdHitSubDetectorNumber[ialgo];
00177 for ( unsigned idet=0; idet<thirdHitSubDetectorNumber[ialgo]; ++idet ) {
00178 thirdHitSubDetectors[ialgo].push_back(thirdSubDets[isub3++]);
00179 }
00180 }
00181 if ( thirdSubDets.size() != check3 )
00182 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00183 << " WARNING : thirdHitSubDetectors does not have the proper size (should be " << check3 << ")"
00184 << std::endl;
00185
00186 originRadius = conf.getParameter<std::vector<double> >("originRadius");
00187 if ( originRadius.size() != seedingAlgo.size() )
00188 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00189 << " WARNING : originRadius does not have the proper size "
00190 << std::endl;
00191
00192 originHalfLength = conf.getParameter<std::vector<double> >("originHalfLength");
00193 if ( originHalfLength.size() != seedingAlgo.size() )
00194 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00195 << " WARNING : originHalfLength does not have the proper size "
00196 << std::endl;
00197
00198 originpTMin = conf.getParameter<std::vector<double> >("originpTMin");
00199 if ( originpTMin.size() != seedingAlgo.size() )
00200 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00201 << " WARNING : originpTMin does not have the proper size "
00202 << std::endl;
00203
00204 primaryVertices = conf.getParameter<std::vector<edm::InputTag> >("primaryVertices");
00205 if ( primaryVertices.size() != seedingAlgo.size() )
00206 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00207 << " WARNING : primaryVertices does not have the proper size "
00208 << std::endl;
00209
00210 zVertexConstraint = conf.getParameter<std::vector<double> >("zVertexConstraint");
00211 if ( zVertexConstraint.size() != seedingAlgo.size() )
00212 throw cms::Exception("FastSimulation/TrajectorySeedProducer ")
00213 << " WARNING : zVertexConstraint does not have the proper size "
00214 << std::endl;
00215 }
00216
00217 }
00218
00219
00220
00221 TrajectorySeedProducer::~TrajectorySeedProducer() {
00222
00223 if(thePropagator) delete thePropagator;
00224
00225
00226 #ifdef FAMOS_DEBUG
00227 std::cout << "TrajectorySeedProducer destructed" << std::endl;
00228 #endif
00229
00230 }
00231
00232 void
00233 TrajectorySeedProducer::beginRun(edm::Run & run, const edm::EventSetup & es) {
00234
00235
00236
00237
00238 edm::ESHandle<MagneticField> magField;
00239 edm::ESHandle<TrackerGeometry> geometry;
00240 edm::ESHandle<MagneticFieldMap> magFieldMap;
00241
00242
00243 es.get<IdealMagneticFieldRecord>().get(magField);
00244 es.get<TrackerDigiGeometryRecord>().get(geometry);
00245 es.get<MagneticFieldMapRecord>().get(magFieldMap);
00246
00247 theMagField = &(*magField);
00248 theGeometry = &(*geometry);
00249 theFieldMap = &(*magFieldMap);
00250
00251 thePropagator = new PropagatorWithMaterial(alongMomentum,0.105,&(*theMagField));
00252
00253 const GlobalPoint g(0.,0.,0.);
00254
00255 }
00256
00257
00258 void
00259 TrajectorySeedProducer::produce(edm::Event& e, const edm::EventSetup& es) {
00260
00261
00262
00263
00264 #ifdef FAMOS_DEBUG
00265 std::cout << "################################################################" << std::endl;
00266 std::cout << " TrajectorySeedProducer produce init " << std::endl;
00267 #endif
00268
00269
00270
00271 unsigned nSimTracks = 0;
00272 unsigned nTracksWithHits = 0;
00273 unsigned nTracksWithPT = 0;
00274 unsigned nTracksWithD0Z0 = 0;
00275
00276 PTrajectoryStateOnDet initialState;
00277
00278
00279 std::vector<TrajectorySeedCollection*>
00280 output(seedingAlgo.size(),static_cast<TrajectorySeedCollection*>(0));
00281 for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) {
00282
00283 output[ialgo] = new TrajectorySeedCollection;
00284 }
00285
00286
00287 edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
00288 e.getByLabel(theBeamSpot,recoBeamSpotHandle);
00289 math::XYZPoint BSPosition_ = recoBeamSpotHandle->position();
00290
00291
00292
00293
00294
00295
00296 x0 = BSPosition_.X();
00297 y0 = BSPosition_.Y();
00298 z0 = BSPosition_.Z();
00299
00300
00301 edm::Handle<edm::SimTrackContainer> theSimTracks;
00302 e.getByLabel("famosSimHits",theSimTracks);
00303
00304 edm::Handle<edm::SimVertexContainer> theSimVtx;
00305 e.getByLabel("famosSimHits",theSimVtx);
00306
00307 #ifdef FAMOS_DEBUG
00308 std::cout << " Step A: SimTracks found " << theSimTracks->size() << std::endl;
00309 #endif
00310
00311
00312 edm::Handle<SiTrackerGSMatchedRecHit2DCollection> theGSRecHits;
00313 e.getByLabel(hitProducer, theGSRecHits);
00314
00315
00316 #ifdef FAMOS_DEBUG
00317 std::cout << " Step B: Full GS RecHits found " << theGSRecHits->size() << std::endl;
00318 #endif
00319 if(theGSRecHits->size() == 0) {
00320 for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) {
00321 std::auto_ptr<TrajectorySeedCollection> p(output[ialgo]);
00322 e.put(p,seedingAlgo[ialgo]);
00323 }
00324 return;
00325 }
00326
00327
00328 vertices = std::vector<const reco::VertexCollection*>
00329 (seedingAlgo.size(),static_cast<const reco::VertexCollection*>(0));
00330 for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) {
00331
00332
00333
00334 edm::Handle<reco::VertexCollection> aHandle;
00335 bool isVertexCollection = e.getByLabel(primaryVertices[ialgo],aHandle);
00336 if (!isVertexCollection ) continue;
00337 vertices[ialgo] = &(*aHandle);
00338 }
00339
00340 #ifdef FAMOS_DEBUG
00341 std::cout << " Step C: Loop over the RecHits, track by track " << std::endl;
00342 #endif
00343
00344
00345 const std::vector<unsigned> theSimTrackIds = theGSRecHits->ids();
00346
00347
00348 for ( unsigned tkId=0; tkId != theSimTrackIds.size(); ++tkId ) {
00349
00350 #ifdef FAMOS_DEBUG
00351 std::cout << "Track number " << tkId << std::endl;
00352 #endif
00353
00354 ++nSimTracks;
00355 unsigned simTrackId = theSimTrackIds[tkId];
00356 const SimTrack& theSimTrack = (*theSimTracks)[simTrackId];
00357 #ifdef FAMOS_DEBUG
00358 std::cout << "Pt = " << std::sqrt(theSimTrack.momentum().Perp2())
00359 << " eta " << theSimTrack.momentum().Eta()
00360 << " pdg ID " << theSimTrack.type()
00361 << std::endl;
00362 #endif
00363
00364
00365 if (selectMuons && abs(theSimTrack.type()) != 13) continue;
00366
00367
00368 int vertexIndex = theSimTrack.vertIndex();
00369 const SimVertex& theSimVertex = (*theSimVtx)[vertexIndex];
00370 #ifdef FAMOS_DEBUG
00371 std::cout << " o SimTrack " << theSimTrack << std::endl;
00372 std::cout << " o SimVertex " << theSimVertex << std::endl;
00373 #endif
00374
00375 BaseParticlePropagator theParticle =
00376 BaseParticlePropagator(
00377 RawParticle(XYZTLorentzVector(theSimTrack.momentum().px(),
00378 theSimTrack.momentum().py(),
00379 theSimTrack.momentum().pz(),
00380 theSimTrack.momentum().e()),
00381 XYZTLorentzVector(theSimVertex.position().x(),
00382 theSimVertex.position().y(),
00383 theSimVertex.position().z(),
00384 theSimVertex.position().t())),
00385 0.,0.,4.);
00386 theParticle.setCharge((*theSimTracks)[simTrackId].charge());
00387
00388 SiTrackerGSMatchedRecHit2DCollection::range theRecHitRange = theGSRecHits->get(simTrackId);
00389 SiTrackerGSMatchedRecHit2DCollection::const_iterator theRecHitRangeIteratorBegin = theRecHitRange.first;
00390 SiTrackerGSMatchedRecHit2DCollection::const_iterator theRecHitRangeIteratorEnd = theRecHitRange.second;
00391 SiTrackerGSMatchedRecHit2DCollection::const_iterator iterRecHit;
00392 SiTrackerGSMatchedRecHit2DCollection::const_iterator iterRecHit1;
00393 SiTrackerGSMatchedRecHit2DCollection::const_iterator iterRecHit2;
00394 SiTrackerGSMatchedRecHit2DCollection::const_iterator iterRecHit3;
00395
00396
00397 unsigned numberOfRecHits = 0;
00398 TrackerRecHit previousHit, currentHit;
00399 for ( iterRecHit = theRecHitRangeIteratorBegin;
00400 iterRecHit != theRecHitRangeIteratorEnd;
00401 ++iterRecHit) {
00402 previousHit = currentHit;
00403 currentHit = TrackerRecHit(&(*iterRecHit),theGeometry);
00404 if ( currentHit.isOnTheSameLayer(previousHit) ) continue;
00405 ++numberOfRecHits;
00406 if ( numberOfRecHits == absMinRecHits ) break;
00407 }
00408
00409
00410 for ( unsigned int ialgo = 0; ialgo < seedingAlgo.size(); ++ialgo ) {
00411
00412 #ifdef FAMOS_DEBUG
00413 std::cout << "Algo " << seedingAlgo[ialgo] << std::endl;
00414 #endif
00415
00416
00417 #ifdef FAMOS_DEBUG
00418 std::cout << "The number of RecHits = " << numberOfRecHits << std::endl;
00419 #endif
00420 if ( numberOfRecHits < minRecHits[ialgo] ) continue;
00421 ++nTracksWithHits;
00422
00423
00424 if ( theSimTrack.momentum().Perp2() < pTMin[ialgo] ) continue;
00425 ++nTracksWithPT;
00426
00427
00428 if ( theParticle.xyImpactParameter(x0,y0) > maxD0[ialgo] ) continue;
00429 if ( fabs( theParticle.zImpactParameter(x0,y0) - z0 ) > maxZ0[ialgo] ) continue;
00430 ++nTracksWithD0Z0;
00431
00432 std::vector<TrackerRecHit >
00433 theSeedHits(numberOfHits[ialgo],
00434 static_cast<TrackerRecHit >(TrackerRecHit()));
00435 TrackerRecHit& theSeedHits0 = theSeedHits[0];
00436 TrackerRecHit& theSeedHits1 = theSeedHits[1];
00437 TrackerRecHit& theSeedHits2 = theSeedHits[2];
00438 bool compatible = false;
00439 for ( iterRecHit1 = theRecHitRangeIteratorBegin; iterRecHit1 != theRecHitRangeIteratorEnd; ++iterRecHit1) {
00440 theSeedHits[0] = TrackerRecHit(&(*iterRecHit1),theGeometry);
00441 #ifdef FAMOS_DEBUG
00442 std::cout << "The first hit position = " << theSeedHits0.globalPosition() << std::endl;
00443 std::cout << "The first hit subDetId = " << theSeedHits0.subDetId() << std::endl;
00444 std::cout << "The first hit layer = " << theSeedHits0.layerNumber() << std::endl;
00445 #endif
00446
00447
00448 bool isInside = true;
00449 if (!selectMuons) {
00450 if (newSyntax)
00451 isInside = true;
00452 else
00453 isInside = theSeedHits0.subDetId() < firstHitSubDetectors[ialgo][0];
00454
00455 if ( isInside ) continue;
00456 }
00457
00458
00459
00460 bool isOndet = true;
00461 if (!selectMuons) {
00462 if (newSyntax)
00463 isOndet = theSeedHits0.isOnRequestedDet(layerList);
00464 else
00465 isOndet = theSeedHits0.isOnRequestedDet(firstHitSubDetectors[ialgo], seedingAlgo[ialgo]);
00466
00467
00468 if ( !isOndet ) continue;
00469 }
00470
00471 #ifdef FAMOS_DEBUG
00472 std::cout << "Apparently the first hit is on the requested detector! " << std::endl;
00473 #endif
00474 for ( iterRecHit2 = iterRecHit1+1; iterRecHit2 != theRecHitRangeIteratorEnd; ++iterRecHit2) {
00475 theSeedHits[1] = TrackerRecHit(&(*iterRecHit2),theGeometry);
00476 #ifdef FAMOS_DEBUG
00477 std::cout << "The second hit position = " << theSeedHits1.globalPosition() << std::endl;
00478 std::cout << "The second hit subDetId = " << theSeedHits1.subDetId() << std::endl;
00479 std::cout << "The second hit layer = " << theSeedHits1.layerNumber() << std::endl;
00480 #endif
00481
00482 if (!selectMuons) {
00483
00484 if (newSyntax)
00485 isInside = true;
00486 else
00487 isInside = theSeedHits1.subDetId() < secondHitSubDetectors[ialgo][0];
00488 if ( isInside ) continue;
00489
00490
00491 if (newSyntax)
00492 isOndet = theSeedHits1.isOnRequestedDet(layerList);
00493 else
00494 isOndet = theSeedHits1.isOnRequestedDet(secondHitSubDetectors[ialgo], seedingAlgo[ialgo]);
00495 if ( !isOndet ) break;
00496 }
00497
00498
00499 if ( theSeedHits1.isOnTheSameLayer(theSeedHits0) ) continue;
00500
00501 #ifdef FAMOS_DEBUG
00502 std::cout << "Apparently the second hit is on the requested detector! " << std::endl;
00503 #endif
00504 GlobalPoint gpos1 = theSeedHits0.globalPosition();
00505 GlobalPoint gpos2 = theSeedHits1.globalPosition();
00506 bool forward = theSeedHits0.isForward();
00507 double error = std::sqrt(theSeedHits0.largerError()+theSeedHits1.largerError());
00508
00509
00510 if(seedingAlgo[0] == "PixelLess" || seedingAlgo[0] == "TobTecLayerPairs"){
00511 compatible = true;
00512
00513 } else {
00514 compatible = compatibleWithBeamAxis(gpos1,gpos2,error,forward,ialgo);
00515 }
00516
00517 #ifdef FAMOS_DEBUG
00518 std::cout << "Algo" << seedingAlgo[0] << "\t Are the two hits compatible with the PV? " << compatible << std::endl;
00519 #endif
00520
00521 if (!selectMuons) {
00522
00523 if ( numberOfHits[ialgo] == 2 ) {
00524
00525 if ( seedingAlgo[0] == "ThirdMixedPairs" ){
00526 compatible = compatible && theSeedHits[0].makesAPairWith3rd(theSeedHits[1]);
00527 } else {
00528 compatible = compatible && theSeedHits[0].makesAPairWith(theSeedHits[1]);
00529
00530
00531
00532
00533
00534
00535
00536 }
00537 }
00538 }
00539
00540
00541 if ( !compatible ) continue;
00542
00543 #ifdef FAMOS_DEBUG
00544 std::cout << "Pair kept! " << std::endl;
00545 #endif
00546
00547
00548 if ( numberOfHits[ialgo] == 2 ) break;
00549
00550 compatible = false;
00551
00552 for ( iterRecHit3 = iterRecHit2+1; iterRecHit3 != theRecHitRangeIteratorEnd; ++iterRecHit3) {
00553 theSeedHits[2] = TrackerRecHit(&(*iterRecHit3),theGeometry);
00554 #ifdef FAMOS_DEBUG
00555 std::cout << "The third hit position = " << theSeedHits2.globalPosition() << std::endl;
00556 std::cout << "The third hit subDetId = " << theSeedHits2.subDetId() << std::endl;
00557 std::cout << "The third hit layer = " << theSeedHits2.layerNumber() << std::endl;
00558 #endif
00559
00560 if (!selectMuons) {
00561
00562 if (newSyntax)
00563 isInside = true;
00564 else
00565 isInside = theSeedHits2.subDetId() < thirdHitSubDetectors[ialgo][0];
00566 if ( isInside ) continue;
00567
00568
00569 if (newSyntax)
00570 isOndet = theSeedHits2.isOnRequestedDet(layerList);
00571 else
00572 isOndet = theSeedHits2.isOnRequestedDet(thirdHitSubDetectors[ialgo], seedingAlgo[ialgo]);
00573
00574 if ( !isOndet ) continue;
00575 }
00576
00577
00578 compatible = !(theSeedHits2.isOnTheSameLayer(theSeedHits1));
00579
00580
00581 if (!selectMuons) compatible = compatible && theSeedHits[0].makesATripletWith(theSeedHits[1],theSeedHits[2]);
00582
00583 #ifdef FAMOS_DEBUG
00584 if ( compatible )
00585 std::cout << "Apparently the third hit is on the requested detector! " << std::endl;
00586 #endif
00587
00588 if ( compatible ) break;
00589
00590 }
00591
00592 if ( compatible ) break;
00593
00594 }
00595
00596 if ( compatible ) break;
00597
00598 }
00599
00600
00601
00602 if ( !compatible ) continue;
00603 #ifdef FAMOS_DEBUG
00604 if ( compatible )
00605 std::cout << "@@@ There is at least a compatible seed" << std::endl;
00606 else
00607 std::cout << "@@@ There is no compatible seed" << std::endl;
00608 #endif
00609
00610
00611 #ifdef FAMOS_DEBUG
00612 std::cout << "Preparing to create the TrajectorySeed" << std::endl;
00613 #endif
00614
00615
00616 edm::OwnVector<TrackingRecHit> recHits;
00617 for ( unsigned ih=0; ih<theSeedHits.size(); ++ih ) {
00618 TrackingRecHit* aTrackingRecHit = theSeedHits[ih].hit()->clone();
00619 recHits.push_back(aTrackingRecHit);
00620 }
00621 #ifdef FAMOS_DEBUG
00622 std::cout << "with " << recHits.size() << " hits." << std::endl;
00623 #endif
00624
00625
00626
00627 GlobalPoint position((*theSimVtx)[vertexIndex].position().x(),
00628 (*theSimVtx)[vertexIndex].position().y(),
00629 (*theSimVtx)[vertexIndex].position().z());
00630
00631
00632 GlobalVector momentum( (*theSimTracks)[simTrackId].momentum().x() ,
00633 (*theSimTracks)[simTrackId].momentum().y() ,
00634 (*theSimTracks)[simTrackId].momentum().z() );
00635
00636 float charge = (*theSimTracks)[simTrackId].charge();
00637
00638 GlobalTrajectoryParameters initialParams(position,momentum,(int)charge,theMagField);
00639
00640 AlgebraicSymMatrix55 errorMatrix= AlgebraicMatrixID();
00641
00642
00643
00644
00645 if(theSeedHits0.subDetId() !=1 || theSeedHits0.subDetId() !=2) errorMatrix = errorMatrix * 0.0000001;
00646
00647
00648
00649 #ifdef FAMOS_DEBUG
00650 std::cout << "TrajectorySeedProducer: SimTrack parameters " << std::endl;
00651 std::cout << "\t\t pT = " << (*theSimTracks)[simTrackId].momentum().Pt() << std::endl;
00652 std::cout << "\t\t eta = " << (*theSimTracks)[simTrackId].momentum().Eta() << std::endl;
00653 std::cout << "\t\t phi = " << (*theSimTracks)[simTrackId].momentum().Phi() << std::endl;
00654 std::cout << "TrajectorySeedProducer: AlgebraicSymMatrix " << errorMatrix << std::endl;
00655 #endif
00656 CurvilinearTrajectoryError initialError(errorMatrix);
00657
00658 FreeTrajectoryState initialFTS(initialParams, initialError);
00659 #ifdef FAMOS_DEBUG
00660 std::cout << "TrajectorySeedProducer: FTS momentum " << initialFTS.momentum() << std::endl;
00661 #endif
00662
00663 const GeomDet* initialLayer = theGeometry->idToDet( recHits.front().geographicalId() );
00664
00665
00666
00667
00668 const TrajectoryStateOnSurface initialTSOS = thePropagator->propagate(initialFTS,initialLayer->surface()) ;
00669 if (!initialTSOS.isValid()) continue;
00670
00671 #ifdef FAMOS_DEBUG
00672 std::cout << "TrajectorySeedProducer: TSOS global momentum " << initialTSOS.globalMomentum() << std::endl;
00673 std::cout << "\t\t\tpT = " << initialTSOS.globalMomentum().perp() << std::endl;
00674 std::cout << "\t\t\teta = " << initialTSOS.globalMomentum().eta() << std::endl;
00675 std::cout << "\t\t\tphi = " << initialTSOS.globalMomentum().phi() << std::endl;
00676 std::cout << "TrajectorySeedProducer: TSOS local momentum " << initialTSOS.localMomentum() << std::endl;
00677 std::cout << "TrajectorySeedProducer: TSOS local error " << initialTSOS.localError().positionError() << std::endl;
00678 std::cout << "TrajectorySeedProducer: TSOS local error matrix " << initialTSOS.localError().matrix() << std::endl;
00679 std::cout << "TrajectorySeedProducer: TSOS surface side " << initialTSOS.surfaceSide() << std::endl;
00680 #endif
00681 stateOnDet(initialTSOS,
00682 recHits.front().geographicalId().rawId(),
00683 initialState);
00684
00685 output[ialgo]->push_back(TrajectorySeed(initialState, recHits, alongMomentum));
00686 #ifdef FAMOS_DEBUG
00687 std::cout << "Trajectory seed created ! " << std::endl;
00688 #endif
00689 break;
00690
00691 }
00692
00693 }
00694
00695 for ( unsigned ialgo=0; ialgo<seedingAlgo.size(); ++ialgo ) {
00696 std::auto_ptr<TrajectorySeedCollection> p(output[ialgo]);
00697 e.put(p,seedingAlgo[ialgo]);
00698 }
00699
00700 }
00701
00702
00703
00704
00705
00706
00707 void
00708 TrajectorySeedProducer::stateOnDet(const TrajectoryStateOnSurface& ts,
00709 unsigned int detid,
00710 PTrajectoryStateOnDet& pts) const
00711 {
00712
00713 const AlgebraicSymMatrix55& m = ts.localError().matrix();
00714
00715 int dim = 5;
00716
00717 float localErrors[15];
00718 int k = 0;
00719 for (int i=0; i<dim; ++i) {
00720 for (int j=0; j<=i; ++j) {
00721 localErrors[k++] = m(i,j);
00722 }
00723 }
00724 int surfaceSide = static_cast<int>(ts.surfaceSide());
00725
00726 pts = PTrajectoryStateOnDet( ts.localParameters(),
00727 localErrors, detid,
00728 surfaceSide);
00729 }
00730
00731 bool
00732 TrajectorySeedProducer::compatibleWithBeamAxis(GlobalPoint& gpos1,
00733 GlobalPoint& gpos2,
00734 double error,
00735 bool forward,
00736 unsigned algo) const {
00737
00738 if ( !seedCleaning ) return true;
00739
00740
00741 XYZTLorentzVector thePos1(gpos1.x(),gpos1.y(),gpos1.z(),0.);
00742 XYZTLorentzVector thePos2(gpos2.x(),gpos2.y(),gpos2.z(),0.);
00743 #ifdef FAMOS_DEBUG
00744 std::cout << "ThePos1 = " << thePos1 << std::endl;
00745 std::cout << "ThePos2 = " << thePos2 << std::endl;
00746 #endif
00747
00748
00749
00750
00751
00752
00753 XYZTLorentzVector theMom2 = (thePos2-thePos1);
00754
00755
00756
00757 ParticlePropagator myPart(theMom2,thePos2,1.,theFieldMap);
00758
00762 bool intersect = myPart.propagateToBeamCylinder(thePos1,originRadius[algo]*1.);
00763 if ( !intersect ) return false;
00764
00765 #ifdef FAMOS_DEBUG
00766 std::cout << "MyPart R = " << myPart.R() << "\t Z = " << myPart.Z()
00767 << "\t pT = " << myPart.Pt() << std::endl;
00768 #endif
00769
00770
00771
00772 if ( myPart.Pt() < originpTMin[algo] ) return false;
00773
00774
00775 if ( fabs(myPart.Z()-z0) > originHalfLength[algo] ) return false;
00776
00777
00778 const reco::VertexCollection* theVertices = vertices[algo];
00779 if (!theVertices) return true;
00780 unsigned nVertices = theVertices->size();
00781 if ( !nVertices || zVertexConstraint[algo] < 0. ) return true;
00782
00783 double R1 = std::sqrt ( (thePos1.X()-x0)*(thePos1.X()-x0)
00784 + (thePos1.Y()-y0)*(thePos1.Y()-y0) );
00785 double R2 = std::sqrt ( (thePos2.X()-x0)*(thePos2.X()-x0)
00786 + (thePos2.Y()-y0)*(thePos2.Y()-y0) );
00787
00788 for ( unsigned iv=0; iv<nVertices; ++iv ) {
00789
00790 double zV = (*theVertices)[iv].z();
00791
00792 double checkRZ1 = forward ?
00793 (thePos1.Z()-zV+zVertexConstraint[algo]) / (thePos2.Z()-zV+zVertexConstraint[algo]) * R2 :
00794 -zVertexConstraint[algo] + R1/R2*(thePos2.Z()-zV+zVertexConstraint[algo]);
00795 double checkRZ2 = forward ?
00796 (thePos1.Z()-zV-zVertexConstraint[algo])/(thePos2.Z()-zV-zVertexConstraint[algo]) * R2 :
00797 +zVertexConstraint[algo] + R1/R2*(thePos2.Z()-zV-zVertexConstraint[algo]);
00798 double checkRZmin = std::min(checkRZ1,checkRZ2)-3.*error;
00799 double checkRZmax = std::max(checkRZ1,checkRZ2)+3.*error;
00800
00801 bool compat = forward ?
00802 checkRZmin < R1 && R1 < checkRZmax :
00803 checkRZmin < thePos1.Z()-zV && thePos1.Z()-zV < checkRZmax;
00804
00805 if ( compat ) return compat;
00806 }
00807
00808 return false;
00809
00810 }
00811