CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/SimTracker/TrackerHitAssociation/src/TrackerHitAssociator.cc

Go to the documentation of this file.
00001 // File: TrackerHitAssociator.cc
00002 
00003 #include <memory>
00004 #include <string>
00005 #include <vector>
00006 
00007 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
00008 //--- for Geometry:
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 //for accumulate
00017 #include <numeric>
00018 #include <iostream>
00019 
00020 using namespace std;
00021 using namespace edm;
00022 
00023 //
00024 // Constructor 
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   // Take by default all tracker SimHits
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   // Step A: Get Inputs
00049   //MAKE THIS PRIVATE MEMBERS 
00050   //  edm::Handle<CrossingFrame<PSimHit> > cf_simhit;
00051   //  std::vector<const CrossingFrame<PSimHit> *> cf_simhitvec;
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  //Loop on PSimHit
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 // Constructor with configurables
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   //if track association there is no need to acces the CrossingFrame
00084   if(!doTrackAssoc_) {
00085     
00086   trackerContainers.clear();
00087   trackerContainers = conf.getParameter<std::vector<std::string> >("ROUList");
00088 
00089     // Step A: Get Inputs
00090     //    edm::Handle<CrossingFrame<PSimHit> > cf_simhit;
00091     //    std::vector<const CrossingFrame<PSimHit> *> cf_simhitvec;
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     //  std::cout << "SIMHITVEC SIZE = " <<  cf_simhitvec.size() << std::endl;
00098     
00099     //    TrackerHits = new MixCollection<PSimHit>(cf_simhitvec);
00100     //std::auto_ptr<MixCollection<PSimHit> > allTrackerHits(TrackerHits);
00101     //   std::auto_ptr<MixCollection<PSimHit> > allTrackerHits(new MixCollection<PSimHit>(cf_simhitvec));
00102     std::auto_ptr<MixCollection<PSimHit> > allTrackerHits(new MixCollection<PSimHit>(cf_simhitvec));
00103     TrackerHits = (*allTrackerHits);
00104     
00105     //Loop on PSimHit
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   //check in case of TTRH
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   //vector with the matched SimHit
00130   std::vector<PSimHit> result; 
00131   //  std::vector<PSimHit> result_old; 
00132   //initialize vectors!
00133   simtrackid.clear();
00134   simhitCFPos.clear();
00135   StripHits = false;
00136   
00137   //get the Detector type of the rechit
00138   DetId detid=  thit.geographicalId();
00139   uint32_t detID = detid.rawId();
00140 
00141   //  cout << "Associator ---> get Detid " << detID << endl;
00142   //check we are in the strip tracker
00143   if(detid.subdetId() == StripSubdetector::TIB ||
00144      detid.subdetId() == StripSubdetector::TOB || 
00145      detid.subdetId() == StripSubdetector::TID ||
00146      detid.subdetId() == StripSubdetector::TEC) 
00147     {
00148       //check if it is a simple SiStripRecHit2D
00149       if(const SiStripRecHit2D * rechit = 
00150          dynamic_cast<const SiStripRecHit2D *>(&thit))
00151         {         
00152           //std::cout << "associate to hit2D" << std::endl;
00153           associateSimpleRecHit(rechit, simtrackid);
00154         }
00155       else if(const SiStripRecHit1D * rechit = 
00156               dynamic_cast<const SiStripRecHit1D *>(&thit)) //check if it is a SiStripRecHit1D
00157         {         
00158           //std::cout << "associate to hit1D" << std::endl;
00159           associateSiStripRecHit1D(rechit,simtrackid);
00160         }
00161       else if(const SiStripMatchedRecHit2D * rechit = 
00162               dynamic_cast<const SiStripMatchedRecHit2D *>(&thit)) //check if it is a matched SiStripMatchedRecHit2D
00163         {         
00164           //std::cout << "associate to matched" << std::endl;
00165           simtrackid = associateMatchedRecHit(rechit);
00166         }
00167       else if(const ProjectedSiStripRecHit2D * rechit = 
00168               dynamic_cast<const ProjectedSiStripRecHit2D *>(&thit)) //check if it is a  ProjectedSiStripRecHit2D
00169         {         
00170           //std::cout << "associate to projectedHit" << std::endl;
00171           simtrackid = associateProjectedRecHit(rechit);
00172           detid = rechit->originalHit().geographicalId();
00173           detID = detid.rawId();
00174         }
00175       else {
00176         //std::cout << "associate to Invalid strip hit??" << std::endl;
00177         //throw cms::Exception("Unknown RecHit Type") << "TrackerHitAssociator failed first casting of " << typeid(thit).name() << " type ";
00178       }
00179 
00180     }
00181   //check we are in the pixel tracker
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           //std::cout << "associate to pixelHit" << std::endl;
00188           associatePixelRecHit(rechit,simtrackid );
00189         }
00190     }
00191   //check if these are GSRecHits (from FastSim)
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   //check if these are GSMatchedRecHits (from FastSim)
00202   if(const SiTrackerGSMatchedRecHit2D * rechit = dynamic_cast<const SiTrackerGSMatchedRecHit2D *>(&thit))
00203     {
00204       simtrackid = associateGSMatchedRecHit(rechit);
00205     }
00206   
00207   //
00208   //Save the SimHits in a vector. for the macthed hits both the rphi and stereo simhits are saved. 
00209   //
00210   
00211   if(StripHits){
00212     //USE THIS FOR STRIPS
00213     //  std::cout << "NEW SIZE =  " << simhitCFPos.size() << std::endl;
00214     
00215     for(size_t i=0; i<simhitCFPos.size(); i++){
00216       //std::cout << "NEW CFPOS " << simhitCFPos[i] << endl;
00217       //std::cout << "NEW LOCALPOS " <<  TrackerHits.getObject(simhitCFPos[i]).localPosition()  << endl;
00218       result.push_back( TrackerHits.getObject(simhitCFPos[i]));
00219     }
00220   }else {
00221     
00222     //now get the SimHit from the trackid
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             //    cout << "Associator ---> ID" << ihit.trackId() << " Simhit x= " << ihit.localPosition().x() 
00238             //         << " y= " <<  ihit.localPosition().y() << " z= " <<  ihit.localPosition().x() << endl; 
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);//iterator to the simhit in the rphi module
00247       std::map<unsigned int, std::vector<PSimHit> >::const_iterator itster = 
00248         SimHitMap.find(detID+1);//iterator to the simhit in the stereo module
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           //===>>>>>change here!!!!
00260           //    for(size_t i=0; i<simtrackid.size();i++){
00261           //  cout << " GluedDet Associator -->  check sihit id's = " << simHitid <<"; compared id's = "<< simtrackid[i] <<endl;
00262           // if(simHitid == simtrackid[i]){ //exclude the geant particles. they all have the same id
00263           for(size_t i=0; i<simtrackid.size();i++){
00264             if(simHitid == simtrackid[i].first && simHiteid == simtrackid[i].second){ 
00265               //          cout << "GluedDet Associator ---> ID" << ihit.trackId() << " Simhit x= " << ihit.localPosition().x() 
00266               //               << " y= " <<  ihit.localPosition().y() << " z= " <<  ihit.localPosition().x() << endl; 
00267               result.push_back(ihit);
00268             }
00269           }
00270         }
00271       }
00272     }
00273   }
00274 
00275 
00276   return result;  
00277 }
00278 
00279 //std::vector<unsigned int> TrackerHitAssociator::associateHitId(const TrackingRecHit & thit) 
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   //check in case of TTRH
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     //vector with the associated Simtkid 
00298     //vector with the CF position of the associated simhits
00299     
00300   //get the Detector type of the rechit
00301     DetId detid=  thit.geographicalId();
00302     //apparently not used
00303     //  uint32_t detID = detid.rawId();
00304     if (const SiTrackerMultiRecHit * rechit = dynamic_cast<const SiTrackerMultiRecHit *>(&thit)){
00305        simtkid=associateMultiRecHitId(rechit);
00306     }
00307     
00308   //cout << "Associator ---> get Detid " << detID << endl;
00309   //check we are in the strip tracker
00310     if(detid.subdetId() == StripSubdetector::TIB ||
00311        detid.subdetId() == StripSubdetector::TOB || 
00312        detid.subdetId() == StripSubdetector::TID ||
00313        detid.subdetId() == StripSubdetector::TEC) 
00314       {
00315         //check if it is a simple SiStripRecHit2D
00316         if(const SiStripRecHit2D * rechit = 
00317            dynamic_cast<const SiStripRecHit2D *>(&thit))
00318           {       
00319             associateSimpleRecHit(rechit, simtkid);
00320           }
00321         //check if it is a matched SiStripMatchedRecHit2D
00322         else  if(const SiStripRecHit1D * rechit = 
00323                  dynamic_cast<const SiStripRecHit1D *>(&thit))
00324           {       
00325             associateSiStripRecHit1D(rechit,simtkid);
00326           }
00327         //check if it is a matched SiStripMatchedRecHit2D
00328         else  if(const SiStripMatchedRecHit2D * rechit = 
00329                  dynamic_cast<const SiStripMatchedRecHit2D *>(&thit))
00330           {       
00331             simtkid = associateMatchedRecHit(rechit);
00332           }
00333         //check if it is a  ProjectedSiStripRecHit2D
00334         else if(const ProjectedSiStripRecHit2D * rechit = 
00335                 dynamic_cast<const ProjectedSiStripRecHit2D *>(&thit))
00336           {       
00337             simtkid = associateProjectedRecHit(rechit);
00338           }
00339         else{
00340           //std::cout << "associate to invalid" << std::endl;
00341           //throw cms::Exception("Unknown RecHit Type") << "TrackerHitAssociator failed second casting of " << typeid(thit).name() << " type ";
00342         }
00343       }
00344     //check we are in the pixel tracker
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     //check if these are GSRecHits (from FastSim)
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 //This function could be templated to avoid to repeat the same code twice??
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 // Caller needs to clear simhit before calling this function
00401 
00402   //initialize vector
00403   theSimtrackid.clear();
00404   //initialize class vector
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 // Caller needs to clear simhit before calling this function
00417 
00418   //initialize class vectors
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   //  std::cout <<"ASSOCIATE SIMPLE RECHIT" << std::endl;           
00431   StripHits =true;        
00432   
00433   DetId detid=  clust->geographicalId();
00434   uint32_t detID = detid.rawId();
00435   
00436   //to store temporary charge information
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()) {  //if it is not empty
00445     //link_detset is a structure, link_detset.data is a std::vector<StripDigiSimLink>
00446     //edm::DetSet<StripDigiSimLink> link_detset = (*stripdigisimlink)[detID];
00447     edm::DetSet<StripDigiSimLink> link_detset = (*isearch);
00448     
00449     if(clust!=0){//the cluster is valid
00450 
00451       //    const edm::Ref<edm::DetSetVector<SiStripCluster>, SiStripCluster, edm::refhelper::FindForDetSetVector<SiStripCluster> > clust=simplerechit->cluster();
00452       
00453       //float chg;
00454       int clusiz = clust->amplitudes().size();
00455       int first  = clust->firstStrip();     
00456       int last   = first + clusiz;
00457       
00458 //       std::cout << "CLUSTERSIZE " << clusiz << " first strip = " << first << " last strip = " << last-1 << std::endl;
00459 //       std::cout << " detID = " << detID << " DETSET size = " << link_detset.data.size() << std::endl;
00460       //use a vector
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         //StripDigiSimLink link = *linkiter;
00465         
00466         if( (int)(linkiter->channel()) >= first  && (int)(linkiter->channel()) < last ){
00467           
00468           //check this digisimlink
00469 //        printf("%s%4d%s%8d%s%3d%s%8.4f\n", "CHANNEL = ", linkiter->channel(), " TrackID = ", linkiter->SimTrackId(),
00470 //               " Process = ", TrackerHits.getObject(linkiter->CFposition()-1).processType(), " fraction = ", linkiter->fraction());
00471           /*
00472             std::cout << "CHECKING CHANNEL  = " << linkiter->channel()   << std::endl;
00473             std::cout << "TrackID  = " << linkiter->SimTrackId()  << std::endl;
00474             std::cout << "Position = " << linkiter->CFposition()  << std::endl;
00475             std::cout << " POS -1 = " << TrackerHits.getObject(linkiter->CFposition()-1).localPosition() << std::endl;
00476             std::cout << " Process = " << TrackerHits.getObject(linkiter->CFposition()-1).processType() << std::endl;
00477             std::cout << " fraction = " << linkiter->fraction() << std::endl;
00478           */
00479           
00480           SimHitIdpr currentId(linkiter->SimTrackId(), linkiter->eventId());
00481           
00482           //create a vector with the list of SimTrack ID's of the tracks that contributed to the RecHit
00483           //write the id only once in the vector
00484           
00485           if(find(idcachev.begin(),idcachev.end(),currentId ) == idcachev.end()){
00486             /*
00487               std::cout << " Adding track id  = " << currentId.first  
00488               << " Event id = " << currentId.second.event() 
00489               << " Bunch Xing = " << currentId.second.bunchCrossing() 
00490               << std::endl;
00491             */
00492             idcachev.push_back(currentId);
00493             simtrackid.push_back(currentId);
00494           }
00495           
00496           //create a vector that contains all the position (in the MixCollection) of the SimHits that contributed to the RecHit
00497           //write position only once
00498           int currentCFPos = linkiter->CFposition()-1;
00499           if(find(CFposcachev.begin(),CFposcachev.end(),currentCFPos ) == CFposcachev.end()){
00500             /*
00501               std::cout << "CHECKING CHANNEL  = " << linkiter->channel()   << std::endl;
00502               std::cout << "\tTrackID  = " << linkiter->SimTrackId()  << "\tCFPos = " << currentCFPos  << std::endl;
00503               std::cout << "\tLocal Pos = " << TrackerHits.getObject(currentCFPos).localPosition() 
00504               << "\tProcess = " << TrackerHits.getObject(currentCFPos).processType() << std::endl;
00505             */
00506             CFposcachev.push_back(currentCFPos);
00507             simhitCFPos.push_back(currentCFPos);
00508             //    simhitassoc.push_back( TrackerHits.getObject(currentCFPos));
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   //associate the two simple hits separately
00533   associateSimpleRecHit(mono,matched_mono );
00534   associateSimpleRecHit(st, matched_st );
00535   
00536   //save in a vector all the simtrack-id's that are common to mono and stereo hits
00537   if(!matched_mono.empty() && !matched_st.empty()){
00538     simtrackid.clear(); //final result vector
00539     //    std::vector<unsigned int> idcachev;
00540     std::vector<SimHitIdpr> idcachev;
00541     //for(vector<unsigned int>::iterator mhit=matched_mono.begin(); mhit != matched_mono.end(); mhit++){
00542     for(vector<SimHitIdpr>::iterator mhit=matched_mono.begin(); mhit != matched_mono.end(); mhit++){
00543       //save only once the ID
00544       if(find(idcachev.begin(), idcachev.end(),(*mhit)) == idcachev.end()) {
00545         idcachev.push_back(*mhit);
00546         //save if the stereoID matched the monoID
00547         if(find(matched_st.begin(), matched_st.end(),(*mhit))!=matched_st.end()){
00548           simtrackid.push_back(*mhit);
00549           //std::cout << "matched case: saved ID " << (*mhit) << std::endl; 
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   //projectedRecHit is a "matched" rechit with only one component
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 //std::vector<unsigned int>  TrackerHitAssociator::associatePixelRecHit(const SiPixelRecHit * pixelrechit)
00574 void  TrackerHitAssociator::associatePixelRecHit(const SiPixelRecHit * pixelrechit, std::vector<SimHitIdpr> & simtrackid)
00575 {
00576   StripHits = false;
00577 
00578   //
00579   // Pixel associator
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()) {  //if it is not empty
00585     //    edm::DetSet<PixelDigiSimLink> link_detset = (*pixeldigisimlink)[detID];
00586     edm::DetSet<PixelDigiSimLink> link_detset = (*isearch);
00587     //    edm::Ref< edm::DetSetVector<SiPixelCluster>, SiPixelCluster> const& cluster = pixelrechit->cluster();
00588     SiPixelRecHit::ClusterRef const& cluster = pixelrechit->cluster();
00589     
00590     //check the reference is valid
00591     
00592     if(!(cluster.isNull())){//if the cluster is valid
00593       
00594       int minPixelRow = (*cluster).minPixelRow();
00595       int maxPixelRow = (*cluster).maxPixelRow();
00596       int minPixelCol = (*cluster).minPixelCol();
00597       int maxPixelCol = (*cluster).maxPixelCol();    
00598       //std::cout << "    Cluster minRow " << minPixelRow << " maxRow " << maxPixelRow << std::endl;
00599       //std::cout << "    Cluster minCol " << minPixelCol << " maxCol " << maxPixelCol << std::endl;
00600       edm::DetSet<PixelDigiSimLink>::const_iterator linkiter = link_detset.data.begin();
00601       int dsl = 0;
00602       //    std::vector<unsigned int> idcachev;
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         //std::cout << "    " << dsl << ") Digi link: row " << pixel_coord.first << " col " << pixel_coord.second << std::endl;      
00608         if(  pixel_coord.first  <= maxPixelRow && 
00609              pixel_coord.first  >= minPixelRow &&
00610              pixel_coord.second <= maxPixelCol &&
00611              pixel_coord.second >= minPixelCol ) {
00612           //std::cout << "      !-> trackid   " << linkiter->SimTrackId() << endl;
00613           //std::cout << "          fraction  " << linkiter->fraction()   << endl;
00614           SimHitIdpr currentId(linkiter->SimTrackId(), linkiter->eventId());
00615           //    if(find(idcachev.begin(),idcachev.end(),linkiter->SimTrackId()) == idcachev.end()){
00616           if(find(idcachev.begin(),idcachev.end(),currentId) == idcachev.end()){
00617             //    simtrackid.push_back(linkiter->SimTrackId());
00618             //idcachev.push_back(linkiter->SimTrackId());
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   //GSRecHit is the FastSimulation RecHit that contains the TrackId already
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   //        std::vector<PSimHit> assimhits;
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   //    std::vector<PSimHit> assimhits = associateHit(**mostpropable);
00656   //    assimhits.insert(assimhits.end(), asstocurrent.begin(), asstocurrent.end());
00657   //std::cout << "Returning " << assimhits.size() << " simhits" << std::endl;
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   //GSRecHit is the FastSimulation RecHit that contains the TrackId already
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