Reference Manual

Go to the documentation of this file.
00001 //
00002 // modified & integrated by Giovanni Abbiendi
00003 // from code by Arun Luthra: UserCode/luthra/MuonTrackSelector/src/MuonTrackSelector.cc
00004 //
00005 #include "SimMuon/MCTruth/plugins/MuonTrackProducer.h"
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00008 #include "DataFormats/TrackReco/interface/Track.h"
00009 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
00010 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
00011 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
00012 #include "DataFormats/MuonReco/interface/MuonSelectors.h"
00013 #include <sstream>
00014 
00015 MuonTrackProducer::MuonTrackProducer(const edm::ParameterSet& parset) :
00016   muonsTag(parset.getParameter< edm::InputTag >("muonsTag")),
00017   inputDTRecSegment4DCollection_(parset.getParameter<edm::InputTag>("inputDTRecSegment4DCollection")),
00018   inputCSCSegmentCollection_(parset.getParameter<edm::InputTag>("inputCSCSegmentCollection")),
00019   selectionTags(parset.getParameter< std::vector<std::string> >("selectionTags")),
00020   trackType(parset.getParameter< std::string >("trackType")),
00021   parset_(parset)
00022 {
00023   edm::LogVerbatim("MuonTrackProducer") << "constructing  MuonTrackProducer" << parset_.dump();
00024   produces<reco::TrackCollection>();
00025   produces<reco::TrackExtraCollection>();
00026   produces<TrackingRecHitCollection>();
00027 }
00028 
00029 MuonTrackProducer::~MuonTrackProducer() {
00030 }
00031 
00032 void MuonTrackProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) 
00033 {
00034   iEvent.getByLabel(muonsTag,muonCollectionH);
00035   iEvent.getByLabel(inputDTRecSegment4DCollection_, dtSegmentCollectionH_);
00036   iEvent.getByLabel(inputCSCSegmentCollection_, cscSegmentCollectionH_);
00037   
00038   std::auto_ptr<reco::TrackCollection> selectedTracks(new reco::TrackCollection);
00039   std::auto_ptr<reco::TrackExtraCollection> selectedTrackExtras( new reco::TrackExtraCollection() );
00040   std::auto_ptr<TrackingRecHitCollection> selectedTrackHits( new TrackingRecHitCollection() );
00041 
00042   reco::TrackRefProd rTracks = iEvent.getRefBeforePut<reco::TrackCollection>();
00043   reco::TrackExtraRefProd rTrackExtras = iEvent.getRefBeforePut<reco::TrackExtraCollection>();
00044   TrackingRecHitRefProd rHits = iEvent.getRefBeforePut<TrackingRecHitCollection>();
00045 
00046   edm::Ref<reco::TrackExtraCollection>::key_type idx = 0;
00047   edm::Ref<reco::TrackExtraCollection>::key_type hidx = 0;
00048 
00049   edm::LogVerbatim("MuonTrackProducer") <<"\nThere are "<< dtSegmentCollectionH_->size()<<" DT segments.";
00050   unsigned int index_dt_segment = 0;
00051   for(DTRecSegment4DCollection::const_iterator segment = dtSegmentCollectionH_->begin();
00052       segment != dtSegmentCollectionH_->end(); ++segment , index_dt_segment++) {
00053     LocalPoint  segmentLocalPosition       = segment->localPosition();
00054     LocalVector segmentLocalDirection      = segment->localDirection();
00055     LocalError  segmentLocalPositionError  = segment->localPositionError();
00056     LocalError  segmentLocalDirectionError = segment->localDirectionError();
00057     DetId geoid = segment->geographicalId();
00058     DTChamberId dtdetid = DTChamberId(geoid);
00059     int wheel = dtdetid.wheel();
00060     int station = dtdetid.station();
00061     int sector = dtdetid.sector();
00062     
00063     float segmentX = segmentLocalPosition.x();
00064     float segmentY = segmentLocalPosition.y();
00065     float segmentdXdZ = segmentLocalDirection.x()/segmentLocalDirection.z();
00066     float segmentdYdZ = segmentLocalDirection.y()/segmentLocalDirection.z();
00067     float segmentXerr = sqrt(segmentLocalPositionError.xx());
00068     float segmentYerr = sqrt(segmentLocalPositionError.yy());
00069     float segmentdXdZerr = sqrt(segmentLocalDirectionError.xx());
00070     float segmentdYdZerr = sqrt(segmentLocalDirectionError.yy());
00071 
00072     edm::LogVerbatim("MuonTrackProducer") 
00073       <<"\nDT segment index :"<<index_dt_segment
00074       <<"\nchamber Wh:"<<wheel<<",St:"<<station<<",Se:"<<sector
00075       <<"\nLocal Position (X,Y)=("<<segmentX<<","<<segmentY<<") +/- ("<<segmentXerr<<","<<segmentYerr<<"), " 
00076       <<"Local Direction (dXdZ,dYdZ)=("<<segmentdXdZ<<","<<segmentdYdZ<<") +/- ("<<segmentdXdZerr<<","<<segmentdYdZerr<<")"; 
00077   }
00078 
00079   edm::LogVerbatim("MuonTrackProducer") <<"\nThere are "<< cscSegmentCollectionH_->size()<<" CSC segments.";
00080   unsigned int index_csc_segment = 0;
00081   for(CSCSegmentCollection::const_iterator segment = cscSegmentCollectionH_->begin();
00082       segment != cscSegmentCollectionH_->end(); ++segment , index_csc_segment++) {
00083     LocalPoint  segmentLocalPosition       = segment->localPosition();
00084     LocalVector segmentLocalDirection      = segment->localDirection();
00085     LocalError  segmentLocalPositionError  = segment->localPositionError();
00086     LocalError  segmentLocalDirectionError = segment->localDirectionError();
00087 
00088     DetId geoid = segment->geographicalId();
00089     CSCDetId cscdetid = CSCDetId(geoid);
00090     int endcap = cscdetid.endcap();
00091     int station = cscdetid.station();
00092     int ring = cscdetid.ring();
00093     int chamber = cscdetid.chamber(); 
00094     
00095     float segmentX = segmentLocalPosition.x();
00096     float segmentY = segmentLocalPosition.y();
00097     float segmentdXdZ = segmentLocalDirection.x()/segmentLocalDirection.z();
00098     float segmentdYdZ = segmentLocalDirection.y()/segmentLocalDirection.z();
00099     float segmentXerr = sqrt(segmentLocalPositionError.xx());
00100     float segmentYerr = sqrt(segmentLocalPositionError.yy());
00101     float segmentdXdZerr = sqrt(segmentLocalDirectionError.xx());
00102     float segmentdYdZerr = sqrt(segmentLocalDirectionError.yy());
00103 
00104     edm::LogVerbatim("MuonTrackProducer") 
00105       <<"\nCSC segment index :"<<index_csc_segment
00106       <<"\nchamber Endcap:"<<endcap<<",St:"<<station<<",Ri:"<<ring<<",Ch:"<<chamber
00107       <<"\nLocal Position (X,Y)=("<<segmentX<<","<<segmentY<<") +/- ("<<segmentXerr<<","<<segmentYerr<<"), " 
00108       <<"Local Direction (dXdZ,dYdZ)=("<<segmentdXdZ<<","<<segmentdYdZ<<") +/- ("<<segmentdXdZerr<<","<<segmentdYdZerr<<")"; 
00109   }
00110 
00111   edm::LogVerbatim("MuonTrackProducer") <<"\nThere are "<< muonCollectionH->size() <<" reco::Muons.";
00112   unsigned int muon_index = 0;
00113   for(reco::MuonCollection::const_iterator muon = muonCollectionH->begin();
00114        muon != muonCollectionH->end(); ++muon, muon_index++) {
00115     edm::LogVerbatim("MuonTrackProducer") <<"\n******* muon index : "<<muon_index;
00116     
00117     std::vector<bool> isGood;    
00118     for(unsigned int index=0; index<selectionTags.size(); ++index) {
00119       isGood.push_back(false);
00120 
00121       muon::SelectionType muonType = muon::selectionTypeFromString(selectionTags[index]);
00122       isGood[index] = muon::isGoodMuon(*muon, muonType);
00123     }
00124 
00125     bool isGoodResult=true;
00126     for(unsigned int index=0; index<isGood.size(); ++index) {
00127       edm::LogVerbatim("MuonTrackProducer") << "selectionTag = "<<selectionTags[index]<< ": "<<isGood[index]<<"\n";
00128       isGoodResult *= isGood[index];
00129     }
00130 
00131     if (isGoodResult) {
00132       // new copy of Track
00133       reco::TrackRef trackref;
00134       if (trackType == "innerTrack") {
00135         if (muon->innerTrack().isNonnull()) trackref = muon->innerTrack();
00136         else continue;
00137       } 
00138       else if (trackType == "outerTrack") {
00139         if (muon->outerTrack().isNonnull()) trackref = muon->outerTrack();
00140         else continue;
00141       } 
00142       else if (trackType == "globalTrack") {
00143         if (muon->globalTrack().isNonnull()) trackref = muon->globalTrack();
00144         else continue;
00145       }
00146       else if (trackType == "innerTrackPlusSegments") {
00147         if (muon->innerTrack().isNonnull()) trackref = muon->innerTrack();
00148         else continue;
00149       }
00150 
00151       const reco::Track* trk = &(*trackref);
00152       // pointer to old track:
00153       reco::Track* newTrk = new reco::Track(*trk);
00154 
00155       newTrk->setExtra( reco::TrackExtraRef( rTrackExtras, idx++ ) );
00156       PropagationDirection seedDir = trk->seedDirection();
00157       // new copy of track Extras
00158       reco::TrackExtra * newExtra = new reco::TrackExtra( trk->outerPosition(), trk->outerMomentum(), 
00159                                         trk->outerOk(), trk->innerPosition(), 
00160                                         trk->innerMomentum(), trk->innerOk(),
00161                                         trk->outerStateCovariance(), trk->outerDetId(),
00162                                         trk->innerStateCovariance(), trk->innerDetId() , seedDir ) ;
00163 
00164       // new copy of the silicon hits; add hit refs to Extra and hits to hit collection
00165       unsigned int index_hit = 0;
00166       
00167       //      edm::LogVerbatim("MuonTrackProducer")<<"\n printing initial hit_pattern";
00168       //      trk->hitPattern().print();
00169         
00170       for (trackingRecHit_iterator iHit = trk->recHitsBegin(); iHit != trk->recHitsEnd(); iHit++) {
00171         TrackingRecHit* hit = (*iHit)->clone();
00172         index_hit++;
00173         selectedTrackHits->push_back( hit );
00174         newExtra->add( TrackingRecHitRef( rHits, hidx++ ) );
00175       }
00176 
00177       if (trackType == "innerTrackPlusSegments") { 
00178         
00179         int wheel, station, sector;
00180         int endcap, /*station, */ ring, chamber;
00181         
00182         edm::LogVerbatim("MuonTrackProducer") <<"Number of chambers: "<<muon->matches().size()
00183                                               <<", arbitrated: "<<muon->numberOfMatches(reco::Muon::SegmentAndTrackArbitration);
00184         unsigned int index_chamber = 0;
00185         
00186         for(std::vector<reco::MuonChamberMatch>::const_iterator chamberMatch = muon->matches().begin();
00187             chamberMatch != muon->matches().end(); ++chamberMatch, index_chamber++) {
00188           std::stringstream chamberStr;
00189           chamberStr <<"\nchamber index: "<<index_chamber; 
00190           
00191           int subdet = chamberMatch->detector();
00192           DetId did = chamberMatch->id;
00193           
00194           if (subdet == MuonSubdetId::DT) {
00195             DTChamberId dtdetid = DTChamberId(did);
00196             wheel = dtdetid.wheel();
00197             station = dtdetid.station();
00198             sector = dtdetid.sector();
00199             chamberStr << ", DT chamber Wh:"<<wheel<<",St:"<<station<<",Se:"<<sector;
00200           } 
00201           else if (subdet == MuonSubdetId::CSC) {
00202             CSCDetId cscdetid = CSCDetId(did);
00203             endcap = cscdetid.endcap();
00204             station = cscdetid.station();
00205             ring = cscdetid.ring();
00206             chamber = cscdetid.chamber();
00207             chamberStr << ", CSC chamber End:"<<endcap<<",St:"<<station<<",Ri:"<<ring<<",Ch:"<<chamber;
00208           }
00209           
00210           chamberStr << ", Number of segments: "<<chamberMatch->segmentMatches.size();
00211           edm::LogVerbatim("MuonTrackProducer") << chamberStr.str();
00212 
00213           unsigned int index_segment = 0;
00214           
00215           for(std::vector<reco::MuonSegmentMatch>::const_iterator segmentMatch = chamberMatch->segmentMatches.begin();
00216               segmentMatch != chamberMatch->segmentMatches.end(); ++segmentMatch, index_segment++) {
00217             
00218             float segmentX = segmentMatch->x;
00219             float segmentY = segmentMatch->y ;
00220             float segmentdXdZ = segmentMatch->dXdZ;
00221             float segmentdYdZ = segmentMatch->dYdZ;
00222             float segmentXerr = segmentMatch->xErr;
00223             float segmentYerr = segmentMatch->yErr;
00224             float segmentdXdZerr = segmentMatch->dXdZErr;
00225             float segmentdYdZerr = segmentMatch->dYdZErr;
00226             
00227             CSCSegmentRef segmentCSC = segmentMatch->cscSegmentRef;
00228             DTRecSegment4DRef segmentDT = segmentMatch->dtSegmentRef;
00229             
00230             bool segment_arbitrated_Ok = (segmentMatch->isMask(reco::MuonSegmentMatch::BestInChamberByDR) && 
00231                                           segmentMatch->isMask(reco::MuonSegmentMatch::BelongsToTrackByDR));
00232             
00233             std::string ARBITRATED(" ***Arbitrated Off*** ");
00234             if (segment_arbitrated_Ok) ARBITRATED = " ***ARBITRATED OK*** ";
00235 
00236             if (subdet == MuonSubdetId::DT) {         
00237               edm::LogVerbatim("MuonTrackProducer")
00238                 <<"\n\t segment index: "<<index_segment << ARBITRATED
00239                 <<"\n\t  Local Position (X,Y)=("<<segmentX<<","<<segmentY<<") +/- ("<<segmentXerr<<","<<segmentYerr<<"), " 
00240                 <<"\n\t  Local Direction (dXdZ,dYdZ)=("<<segmentdXdZ<<","<<segmentdYdZ<<") +/- ("<<segmentdXdZerr<<","<<segmentdYdZerr<<")"; 
00241               
00242               if (!segment_arbitrated_Ok) continue;
00243               
00244               if (segmentDT.get() != 0) {
00245                 const DTRecSegment4D* segment = segmentDT.get();
00246                 
00247                 edm::LogVerbatim("MuonTrackProducer")<<"\t ===> MATCHING with DT segment with index = "<<segmentDT.key();
00248                 
00249                 if(segment->hasPhi()) {
00250                   const DTChamberRecSegment2D* phiSeg = segment->phiSegment();
00251                   std::vector<const TrackingRecHit*> phiHits = phiSeg->recHits();
00252                   for(std::vector<const TrackingRecHit*>::const_iterator ihit = phiHits.begin();
00253                       ihit != phiHits.end(); ++ihit) {
00254                     TrackingRecHit* seghit = (*ihit)->clone();
00255                     newTrk->setHitPattern( *seghit, index_hit);
00256                     //              edm::LogVerbatim("MuonTrackProducer")<<"hit pattern for position "<<index_hit<<" set to:";
00257                     //              newTrk->hitPattern().printHitPattern(index_hit, std::cout);
00258                     index_hit++;
00259                     selectedTrackHits->push_back( seghit );
00260                     newExtra->add( TrackingRecHitRef( rHits, hidx ++ ) );
00261                   }
00262                 }
00263                 
00264                 if(segment->hasZed()) {
00265                   const DTSLRecSegment2D* zSeg = (*segment).zSegment();
00266                   std::vector<const TrackingRecHit*> zedHits = zSeg->recHits();
00267                   for(std::vector<const TrackingRecHit*>::const_iterator ihit = zedHits.begin();
00268                       ihit != zedHits.end(); ++ihit) {
00269                     TrackingRecHit* seghit = (*ihit)->clone();
00270                     newTrk->setHitPattern( *seghit, index_hit);
00271                     //              edm::LogVerbatim("MuonTrackProducer")<<"hit pattern for position "<<index_hit<<" set to:";
00272                     //              newTrk->hitPattern().printHitPattern(index_hit, std::cout);
00273                     index_hit++;
00274                     selectedTrackHits->push_back( seghit );
00275                     newExtra->add( TrackingRecHitRef( rHits, hidx ++ ) );
00276                   }
00277                 }
00278               } else edm::LogWarning("MuonTrackProducer")<<"\n***WARNING: UNMATCHED DT segment ! \n";
00279             } // if (subdet == MuonSubdetId::DT)
00280 
00281             else if (subdet == MuonSubdetId::CSC) {
00282               edm::LogVerbatim("MuonTrackProducer")
00283                 <<"\n\t segment index: "<<index_segment << ARBITRATED
00284                 <<"\n\t  Local Position (X,Y)=("<<segmentX<<","<<segmentY<<") +/- ("<<segmentXerr<<","<<segmentYerr<<"), " 
00285                 <<"\n\t  Local Direction (dXdZ,dYdZ)=("<<segmentdXdZ<<","<<segmentdYdZ<<") +/- ("<<segmentdXdZerr<<","<<segmentdYdZerr<<")"; 
00286               
00287               if (!segment_arbitrated_Ok) continue;
00288               
00289               if (segmentCSC.get() != 0) {
00290                 const CSCSegment* segment = segmentCSC.get();
00291                 
00292                 edm::LogVerbatim("MuonTrackProducer")<<"\t ===> MATCHING with CSC segment with index = "<<segmentCSC.key();
00293                 
00294                 std::vector<const TrackingRecHit*> hits = segment->recHits();
00295                 for(std::vector<const TrackingRecHit*>::const_iterator ihit = hits.begin();
00296                     ihit != hits.end(); ++ihit) {
00297                   TrackingRecHit* seghit = (*ihit)->clone();
00298                   newTrk->setHitPattern( *seghit, index_hit);
00299                   //                edm::LogVerbatim("MuonTrackProducer")<<"hit pattern for position "<<index_hit<<" set to:";
00300                   //                newTrk->hitPattern().printHitPattern(index_hit, std::cout);
00301                   index_hit++;
00302                   selectedTrackHits->push_back( seghit );
00303                   newExtra->add( TrackingRecHitRef( rHits, hidx ++ ) );           
00304                 }
00305               } else edm::LogWarning("MuonTrackProducer")<<"\n***WARNING: UNMATCHED CSC segment ! \n";
00306             }  //  else if (subdet == MuonSubdetId::CSC)
00307 
00308           } // loop on vector<MuonSegmentMatch>   
00309         } // loop on vector<MuonChamberMatch>   
00310       } // if (trackType == "innerTrackPlusSegments")
00311       
00312       //      edm::LogVerbatim("MuonTrackProducer")<<"\n printing final hit_pattern";
00313       //      newTrk->hitPattern().print();
00314       
00315       selectedTracks->push_back( *newTrk );
00316       selectedTrackExtras->push_back( *newExtra );
00317 
00318     } // if (isGoodResult)
00319   }  // loop on reco::MuonCollection
00320   
00321   iEvent.put(selectedTracks);
00322   iEvent.put(selectedTrackExtras);
00323   iEvent.put(selectedTrackHits);
00324 }
/data/refman/pasoursint/CMSSW_5_3_3/src/SimMuon/MCTruth/plugins/MuonTrackProducer.cc
