00001
00002
00003 #include <memory>
00004 #include <string>
00005 #include <vector>
00006
00007 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
00008
00009 #include "DataFormats/Common/interface/Ref.h"
00010 #include "DataFormats/DetId/interface/DetId.h"
00011 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00012 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
00013
00014 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
00015
00016
00017 #include <numeric>
00018 #include <iostream>
00019
00020 using namespace std;
00021 using namespace edm;
00022
00023
00024
00025
00026 TrackerHitAssociator::TrackerHitAssociator(const edm::Event& e) :
00027 myEvent_(e),
00028 doPixel_( true ),
00029 doStrip_( true ),
00030 doTrackAssoc_( false ) {
00031 trackerContainers.clear();
00032
00033
00034
00035 trackerContainers.push_back("g4SimHitsTrackerHitsTIBLowTof");
00036 trackerContainers.push_back("g4SimHitsTrackerHitsTIBHighTof");
00037 trackerContainers.push_back("g4SimHitsTrackerHitsTIDLowTof");
00038 trackerContainers.push_back("g4SimHitsTrackerHitsTIDHighTof");
00039 trackerContainers.push_back("g4SimHitsTrackerHitsTOBLowTof");
00040 trackerContainers.push_back("g4SimHitsTrackerHitsTOBHighTof");
00041 trackerContainers.push_back("g4SimHitsTrackerHitsTECLowTof");
00042 trackerContainers.push_back("g4SimHitsTrackerHitsTECHighTof");
00043 trackerContainers.push_back("g4SimHitsTrackerHitsPixelBarrelLowTof");
00044 trackerContainers.push_back("g4SimHitsTrackerHitsPixelBarrelHighTof");
00045 trackerContainers.push_back("g4SimHitsTrackerHitsPixelEndcapLowTof");
00046 trackerContainers.push_back("g4SimHitsTrackerHitsPixelEndcapHighTof");
00047
00048
00049
00050
00051
00052
00053 for(uint32_t i = 0; i< trackerContainers.size();i++){
00054 e.getByLabel("mix",trackerContainers[i],cf_simhit);
00055 cf_simhitvec.push_back(cf_simhit.product());
00056 }
00057
00058 std::auto_ptr<MixCollection<PSimHit> > allTrackerHits(new MixCollection<PSimHit>(cf_simhitvec));
00059 TrackerHits = (*allTrackerHits);
00060
00061
00062 SimHitMap.clear();
00063
00064 MixCollection<PSimHit>::iterator isim;
00065 for (isim=allTrackerHits->begin(); isim!= allTrackerHits->end();isim++) {
00066 SimHitMap[(*isim).detUnitId()].push_back((*isim));
00067 }
00068
00069 if(doStrip_) e.getByLabel("simSiStripDigis", stripdigisimlink);
00070 if(doPixel_) e.getByLabel("simSiPixelDigis", pixeldigisimlink);
00071
00072 }
00073
00074
00075
00076
00077 TrackerHitAssociator::TrackerHitAssociator(const edm::Event& e, const edm::ParameterSet& conf) :
00078 myEvent_(e),
00079 doPixel_( conf.getParameter<bool>("associatePixel") ),
00080 doStrip_( conf.getParameter<bool>("associateStrip") ),
00081 doTrackAssoc_( conf.getParameter<bool>("associateRecoTracks") ){
00082
00083
00084 if(!doTrackAssoc_) {
00085
00086 trackerContainers.clear();
00087 trackerContainers = conf.getParameter<std::vector<std::string> >("ROUList");
00088
00089
00090
00091
00092 for(uint32_t i = 0; i< trackerContainers.size();i++){
00093 e.getByLabel("mix",trackerContainers[i],cf_simhit);
00094 cf_simhitvec.push_back(cf_simhit.product());
00095 }
00096
00097
00098
00099
00100
00101
00102 std::auto_ptr<MixCollection<PSimHit> > allTrackerHits(new MixCollection<PSimHit>(cf_simhitvec));
00103 TrackerHits = (*allTrackerHits);
00104
00105
00106 SimHitMap.clear();
00107
00108 MixCollection<PSimHit>::iterator isim;
00109 for (isim=allTrackerHits->begin(); isim!= allTrackerHits->end();isim++) {
00110 SimHitMap[(*isim).detUnitId()].push_back((*isim));
00111 }
00112
00113 }
00114
00115 if(doStrip_) e.getByLabel("simSiStripDigis", stripdigisimlink);
00116 if(doPixel_) e.getByLabel("simSiPixelDigis", pixeldigisimlink);
00117
00118 }
00119
00120 std::vector<PSimHit> TrackerHitAssociator::associateHit(const TrackingRecHit & thit)
00121 {
00122
00123
00124 if(const TransientTrackingRecHit * ttrh = dynamic_cast<const TransientTrackingRecHit *>(&thit)) {
00125 std::cout << "calling associateHit for TransientTRH" << std::endl;
00126 return associateHit(*ttrh->hit());
00127 }
00128
00129
00130 std::vector<PSimHit> result;
00131
00132
00133 simtrackid.clear();
00134 simhitCFPos.clear();
00135 StripHits = false;
00136
00137
00138 DetId detid= thit.geographicalId();
00139 uint32_t detID = detid.rawId();
00140
00141
00142
00143 if(detid.subdetId() == StripSubdetector::TIB ||
00144 detid.subdetId() == StripSubdetector::TOB ||
00145 detid.subdetId() == StripSubdetector::TID ||
00146 detid.subdetId() == StripSubdetector::TEC)
00147 {
00148
00149 if(const SiStripRecHit2D * rechit =
00150 dynamic_cast<const SiStripRecHit2D *>(&thit))
00151 {
00152
00153 associateSimpleRecHit(rechit, simtrackid);
00154 }
00155 else if(const SiStripRecHit1D * rechit =
00156 dynamic_cast<const SiStripRecHit1D *>(&thit))
00157 {
00158
00159 associateSiStripRecHit1D(rechit,simtrackid);
00160 }
00161 else if(const SiStripMatchedRecHit2D * rechit =
00162 dynamic_cast<const SiStripMatchedRecHit2D *>(&thit))
00163 {
00164
00165 simtrackid = associateMatchedRecHit(rechit);
00166 }
00167 else if(const ProjectedSiStripRecHit2D * rechit =
00168 dynamic_cast<const ProjectedSiStripRecHit2D *>(&thit))
00169 {
00170
00171 simtrackid = associateProjectedRecHit(rechit);
00172 detid = rechit->originalHit().geographicalId();
00173 detID = detid.rawId();
00174 }
00175 else {
00176
00177
00178 }
00179
00180 }
00181
00182 if( (unsigned int)(detid.subdetId()) == PixelSubdetector::PixelBarrel ||
00183 (unsigned int)(detid.subdetId()) == PixelSubdetector::PixelEndcap)
00184 {
00185 if(const SiPixelRecHit * rechit = dynamic_cast<const SiPixelRecHit *>(&thit))
00186 {
00187
00188 associatePixelRecHit(rechit,simtrackid );
00189 }
00190 }
00191
00192
00193 if(const SiTrackerGSRecHit2D * rechit = dynamic_cast<const SiTrackerGSRecHit2D *>(&thit))
00194 {
00195 simtrackid = associateGSRecHit(rechit);
00196 }
00197 if (const SiTrackerMultiRecHit * rechit = dynamic_cast<const SiTrackerMultiRecHit *>(&thit)){
00198 return associateMultiRecHit(rechit);
00199 }
00200
00201
00202 if(const SiTrackerGSMatchedRecHit2D * rechit = dynamic_cast<const SiTrackerGSMatchedRecHit2D *>(&thit))
00203 {
00204 simtrackid = associateGSMatchedRecHit(rechit);
00205 }
00206
00207
00208
00209
00210
00211 if(StripHits){
00212
00213
00214
00215 for(size_t i=0; i<simhitCFPos.size(); i++){
00216
00217
00218 result.push_back( TrackerHits.getObject(simhitCFPos[i]));
00219 }
00220 }else {
00221
00222
00223 vector<PSimHit> simHit;
00224 std::map<unsigned int, std::vector<PSimHit> >::const_iterator it = SimHitMap.find(detID);
00225 simHit.clear();
00226 if (it!= SimHitMap.end()){
00227 simHit = it->second;
00228 vector<PSimHit>::const_iterator simHitIter = simHit.begin();
00229 vector<PSimHit>::const_iterator simHitIterEnd = simHit.end();
00230 for (;simHitIter != simHitIterEnd; ++simHitIter) {
00231 const PSimHit ihit = *simHitIter;
00232 unsigned int simHitid = ihit.trackId();
00233 EncodedEventId simHiteid = ihit.eventId();
00234
00235 for(size_t i=0; i<simtrackid.size();i++){
00236 if(simHitid == simtrackid[i].first && simHiteid == simtrackid[i].second){
00237
00238
00239 result.push_back(ihit);
00240 }
00241 }
00242 }
00243 }else{
00245 std::map<unsigned int, std::vector<PSimHit> >::const_iterator itrphi =
00246 SimHitMap.find(detID+2);
00247 std::map<unsigned int, std::vector<PSimHit> >::const_iterator itster =
00248 SimHitMap.find(detID+1);
00249 if (itrphi!= SimHitMap.end()&&itster!=SimHitMap.end()){
00250 simHit = itrphi->second;
00251 simHit.insert(simHit.end(),(itster->second).begin(),(itster->second).end());
00252 vector<PSimHit>::const_iterator simHitIter = simHit.begin();
00253 vector<PSimHit>::const_iterator simHitIterEnd = simHit.end();
00254 for (;simHitIter != simHitIterEnd; ++simHitIter) {
00255 const PSimHit ihit = *simHitIter;
00256 unsigned int simHitid = ihit.trackId();
00257 EncodedEventId simHiteid = ihit.eventId();
00258
00259
00260
00261
00262
00263 for(size_t i=0; i<simtrackid.size();i++){
00264 if(simHitid == simtrackid[i].first && simHiteid == simtrackid[i].second){
00265
00266
00267 result.push_back(ihit);
00268 }
00269 }
00270 }
00271 }
00272 }
00273 }
00274
00275
00276 return result;
00277 }
00278
00279
00280 std::vector< SimHitIdpr > TrackerHitAssociator::associateHitId(const TrackingRecHit & thit)
00281 {
00282 std::vector< SimHitIdpr > simhitid;
00283 associateHitId(thit, simhitid);
00284 return simhitid;
00285 }
00286
00287 void TrackerHitAssociator::associateHitId(const TrackingRecHit & thit, std::vector< SimHitIdpr > & simtkid)
00288 {
00289
00290
00291 if(const TransientTrackingRecHit * ttrh = dynamic_cast<const TransientTrackingRecHit *>(&thit)) {
00292 associateHitId(*ttrh->hit(), simtkid);
00293 }
00294 else{
00295 simtkid.clear();
00296 simhitCFPos.clear();
00297
00298
00299
00300
00301 DetId detid= thit.geographicalId();
00302
00303
00304 if (const SiTrackerMultiRecHit * rechit = dynamic_cast<const SiTrackerMultiRecHit *>(&thit)){
00305 simtkid=associateMultiRecHitId(rechit);
00306 }
00307
00308
00309
00310 if(detid.subdetId() == StripSubdetector::TIB ||
00311 detid.subdetId() == StripSubdetector::TOB ||
00312 detid.subdetId() == StripSubdetector::TID ||
00313 detid.subdetId() == StripSubdetector::TEC)
00314 {
00315
00316 if(const SiStripRecHit2D * rechit =
00317 dynamic_cast<const SiStripRecHit2D *>(&thit))
00318 {
00319 associateSimpleRecHit(rechit, simtkid);
00320 }
00321
00322 else if(const SiStripRecHit1D * rechit =
00323 dynamic_cast<const SiStripRecHit1D *>(&thit))
00324 {
00325 associateSiStripRecHit1D(rechit,simtkid);
00326 }
00327
00328 else if(const SiStripMatchedRecHit2D * rechit =
00329 dynamic_cast<const SiStripMatchedRecHit2D *>(&thit))
00330 {
00331 simtkid = associateMatchedRecHit(rechit);
00332 }
00333
00334 else if(const ProjectedSiStripRecHit2D * rechit =
00335 dynamic_cast<const ProjectedSiStripRecHit2D *>(&thit))
00336 {
00337 simtkid = associateProjectedRecHit(rechit);
00338 }
00339 else{
00340
00341
00342 }
00343 }
00344
00345 else if( (unsigned int)(detid.subdetId()) == PixelSubdetector::PixelBarrel ||
00346 (unsigned int)(detid.subdetId()) == PixelSubdetector::PixelEndcap)
00347 {
00348 if(const SiPixelRecHit * rechit = dynamic_cast<const SiPixelRecHit *>(&thit))
00349 {
00350 associatePixelRecHit(rechit,simtkid );
00351 }
00352 }
00353
00354 if(const SiTrackerGSRecHit2D * rechit = dynamic_cast<const SiTrackerGSRecHit2D *>(&thit))
00355 {
00356 simtkid = associateGSRecHit(rechit);
00357 }
00358 if(const SiTrackerGSMatchedRecHit2D * rechit = dynamic_cast<const SiTrackerGSMatchedRecHit2D *>(&thit))
00359 {
00360 simtkid = associateGSMatchedRecHit(rechit);
00361 }
00362
00363 }
00364 }
00365
00366
00367 void TrackerHitAssociator::associateSimpleRecHit(const SiStripRecHit2D * simplerechit, std::vector<SimHitIdpr>& simtrackid)
00368 {
00369 const SiStripCluster* clust = 0;
00370 if(simplerechit->cluster().isNonnull())
00371 {
00372 clust=&(*simplerechit->cluster());
00373 }
00374 else if(simplerechit->cluster_regional().isNonnull())
00375 {
00376 clust=&(*simplerechit->cluster_regional());
00377 }
00378
00379 associateSimpleRecHitCluster(clust,simtrackid);
00380 }
00381
00382
00383
00384 void TrackerHitAssociator::associateSiStripRecHit1D(const SiStripRecHit1D * simplerechit, std::vector<SimHitIdpr>& simtrackid)
00385 {
00386 const SiStripCluster* clust = 0;
00387 if(simplerechit->cluster().isNonnull())
00388 {
00389 clust=&(*simplerechit->cluster());
00390 }
00391 else if(simplerechit->cluster_regional().isNonnull())
00392 {
00393 clust=&(*simplerechit->cluster_regional());
00394 }
00395 associateSimpleRecHitCluster(clust,simtrackid);
00396 }
00397
00398 void TrackerHitAssociator::associateSimpleRecHitCluster(const SiStripCluster* clust, std::vector<SimHitIdpr>& theSimtrackid, std::vector<PSimHit>& simhit)
00399 {
00400
00401
00402
00403 theSimtrackid.clear();
00404
00405 simhitCFPos.clear();
00406
00407 associateSimpleRecHitCluster(clust, theSimtrackid);
00408
00409 for(size_t i=0; i<simhitCFPos.size(); i++){
00410 simhit.push_back(TrackerHits.getObject(simhitCFPos[i]));
00411 }
00412 }
00413
00414 void TrackerHitAssociator::associateSimpleRecHitCluster(const SiStripCluster* clust, std::vector<PSimHit>& simhit)
00415 {
00416
00417
00418
00419 simtrackid.clear();
00420 simhitCFPos.clear();
00421
00422 associateSimpleRecHitCluster(clust, simtrackid);
00423
00424 for(size_t i=0; i<simhitCFPos.size(); i++){
00425 simhit.push_back(TrackerHits.getObject(simhitCFPos[i]));
00426 }
00427 }
00428
00429 void TrackerHitAssociator::associateSimpleRecHitCluster(const SiStripCluster* clust, std::vector<SimHitIdpr>& simtrackid){
00430
00431 StripHits =true;
00432
00433 DetId detid= clust->geographicalId();
00434 uint32_t detID = detid.rawId();
00435
00436
00437 std::vector<SimHitIdpr> cache_simtrackid;
00438 cache_simtrackid.clear();
00439
00440 std::map<SimHitIdpr, vector<float> > temp_simtrackid;
00441 temp_simtrackid.clear();
00442
00443 edm::DetSetVector<StripDigiSimLink>::const_iterator isearch = stripdigisimlink->find(detID);
00444 if(isearch != stripdigisimlink->end()) {
00445
00446
00447 edm::DetSet<StripDigiSimLink> link_detset = (*isearch);
00448
00449 if(clust!=0){
00450
00451
00452
00453
00454 int clusiz = clust->amplitudes().size();
00455 int first = clust->firstStrip();
00456 int last = first + clusiz;
00457
00458
00459
00460
00461 std::vector<SimHitIdpr> idcachev;
00462 std::vector<int> CFposcachev;
00463 for(edm::DetSet<StripDigiSimLink>::const_iterator linkiter = link_detset.data.begin(); linkiter != link_detset.data.end(); linkiter++){
00464
00465
00466 if( (int)(linkiter->channel()) >= first && (int)(linkiter->channel()) < last ){
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480 SimHitIdpr currentId(linkiter->SimTrackId(), linkiter->eventId());
00481
00482
00483
00484
00485 if(find(idcachev.begin(),idcachev.end(),currentId ) == idcachev.end()){
00486
00487
00488
00489
00490
00491
00492 idcachev.push_back(currentId);
00493 simtrackid.push_back(currentId);
00494 }
00495
00496
00497
00498 int currentCFPos = linkiter->CFposition()-1;
00499 if(find(CFposcachev.begin(),CFposcachev.end(),currentCFPos ) == CFposcachev.end()){
00500
00501
00502
00503
00504
00505
00506 CFposcachev.push_back(currentCFPos);
00507 simhitCFPos.push_back(currentCFPos);
00508
00509 }
00510 }
00511 }
00512 }
00513 else {
00514 edm::LogError("TrackerHitAssociator")<<"no cluster reference attached";
00515 }
00516 }
00517 }
00518
00519 std::vector<SimHitIdpr> TrackerHitAssociator::associateMatchedRecHit(const SiStripMatchedRecHit2D * matchedrechit)
00520 {
00521
00522 StripHits = true;
00523
00524
00525 vector<SimHitIdpr> matched_mono;
00526 vector<SimHitIdpr> matched_st;
00527 matched_mono.clear();
00528 matched_st.clear();
00529
00530 const SiStripRecHit2D *mono = matchedrechit->monoHit();
00531 const SiStripRecHit2D *st = matchedrechit->stereoHit();
00532
00533 associateSimpleRecHit(mono,matched_mono );
00534 associateSimpleRecHit(st, matched_st );
00535
00536
00537 if(!matched_mono.empty() && !matched_st.empty()){
00538 simtrackid.clear();
00539
00540 std::vector<SimHitIdpr> idcachev;
00541
00542 for(vector<SimHitIdpr>::iterator mhit=matched_mono.begin(); mhit != matched_mono.end(); mhit++){
00543
00544 if(find(idcachev.begin(), idcachev.end(),(*mhit)) == idcachev.end()) {
00545 idcachev.push_back(*mhit);
00546
00547 if(find(matched_st.begin(), matched_st.end(),(*mhit))!=matched_st.end()){
00548 simtrackid.push_back(*mhit);
00549
00550 }
00551 }
00552 }
00553 }
00554 return simtrackid;
00555 }
00556
00557
00558 std::vector<SimHitIdpr> TrackerHitAssociator::associateProjectedRecHit(const ProjectedSiStripRecHit2D * projectedrechit)
00559 {
00560 StripHits = true;
00561
00562
00563
00564
00565 vector<SimHitIdpr> matched_mono;
00566 matched_mono.clear();
00567
00568 const SiStripRecHit2D mono = projectedrechit->originalHit();
00569 associateSimpleRecHit(&mono, matched_mono);
00570 return matched_mono;
00571 }
00572
00573
00574 void TrackerHitAssociator::associatePixelRecHit(const SiPixelRecHit * pixelrechit, std::vector<SimHitIdpr> & simtrackid)
00575 {
00576 StripHits = false;
00577
00578
00579
00580
00581 DetId detid= pixelrechit->geographicalId();
00582 uint32_t detID = detid.rawId();
00583 edm::DetSetVector<PixelDigiSimLink>::const_iterator isearch = pixeldigisimlink->find(detID);
00584 if(isearch != pixeldigisimlink->end()) {
00585
00586 edm::DetSet<PixelDigiSimLink> link_detset = (*isearch);
00587
00588 SiPixelRecHit::ClusterRef const& cluster = pixelrechit->cluster();
00589
00590
00591
00592 if(!(cluster.isNull())){
00593
00594 int minPixelRow = (*cluster).minPixelRow();
00595 int maxPixelRow = (*cluster).maxPixelRow();
00596 int minPixelCol = (*cluster).minPixelCol();
00597 int maxPixelCol = (*cluster).maxPixelCol();
00598
00599
00600 edm::DetSet<PixelDigiSimLink>::const_iterator linkiter = link_detset.data.begin();
00601 int dsl = 0;
00602
00603 std::vector<SimHitIdpr> idcachev;
00604 for( ; linkiter != link_detset.data.end(); linkiter++) {
00605 dsl++;
00606 std::pair<int,int> pixel_coord = PixelDigi::channelToPixel(linkiter->channel());
00607
00608 if( pixel_coord.first <= maxPixelRow &&
00609 pixel_coord.first >= minPixelRow &&
00610 pixel_coord.second <= maxPixelCol &&
00611 pixel_coord.second >= minPixelCol ) {
00612
00613
00614 SimHitIdpr currentId(linkiter->SimTrackId(), linkiter->eventId());
00615
00616 if(find(idcachev.begin(),idcachev.end(),currentId) == idcachev.end()){
00617
00618
00619 simtrackid.push_back(currentId);
00620 idcachev.push_back(currentId);
00621 }
00622 }
00623 }
00624 }
00625 else{
00626 edm::LogError("TrackerHitAssociator")<<"no Pixel cluster reference attached";
00627
00628 }
00629 }
00630
00631
00632 }
00633
00634 std::vector<SimHitIdpr> TrackerHitAssociator::associateGSRecHit(const SiTrackerGSRecHit2D * gsrechit)
00635 {
00636 StripHits = false;
00637
00638
00639 vector<SimHitIdpr> simtrackid;
00640 simtrackid.clear();
00641 SimHitIdpr currentId(gsrechit->simtrackId(), EncodedEventId(gsrechit->eeId()));
00642 simtrackid.push_back(currentId);
00643 return simtrackid;
00644 }
00645
00646 std::vector<PSimHit> TrackerHitAssociator::associateMultiRecHit(const SiTrackerMultiRecHit * multirechit){
00647 std::vector<const TrackingRecHit*> componenthits = multirechit->recHits();
00648
00649 int size=multirechit->weights().size(), idmostprobable=0;
00650
00651 for (int i=0; i<size; i++){
00652 if(multirechit->weight(i)>multirechit->weight(idmostprobable)) idmostprobable=i;
00653 }
00654
00655
00656
00657
00658 return associateHit(*componenthits[idmostprobable]);
00659 }
00660
00661 std::vector<SimHitIdpr> TrackerHitAssociator::associateMultiRecHitId(const SiTrackerMultiRecHit * multirechit){
00662 std::vector<const TrackingRecHit*> componenthits = multirechit->recHits();
00663 int size=multirechit->weights().size(), idmostprobable=0;
00664
00665 for (int i=0; i<size; i++){
00666 if(multirechit->weight(i)>multirechit->weight(idmostprobable)) idmostprobable=i;
00667 }
00668
00669 return associateHitId(*componenthits[idmostprobable]);
00670 }
00671
00672 std::vector<SimHitIdpr> TrackerHitAssociator::associateGSMatchedRecHit(const SiTrackerGSMatchedRecHit2D * gsmrechit)
00673 {
00674 StripHits = false;
00675
00676
00677 vector<SimHitIdpr> simtrackid;
00678 simtrackid.clear();
00679 SimHitIdpr currentId(gsmrechit->simtrackId(), EncodedEventId(gsmrechit->eeId()));
00680 simtrackid.push_back(currentId);
00681 return simtrackid;
00682 }
00683